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ABSTRACT 


An analytical and computerized study of the steady state and 
transient response of a phosphoric acid fuel cell (PAFC) system was 
completed. Parametric studies and sensitivity analyses of the PAFC 
system's operation were accomplished. 

Four non-linear dynamic models of the fuel cell staci;, reformer, 
shift converters, and heat exchangers were developed based on 
nonhomosenaous non-linear partial differential equations, which include 
the material, component, energy balance, and electrochemical kinetic 
features. Due to a lack of experimental data for the dynamic response 
of the coiT.j%nents, only the steady state results were compared with data 
from other sources, indicating reasonably good agreement. 

A steady state simulation of the entire system was developed using 
nonlinear ordinary differential equations. The finite difference method 
and trial-and-error procedures were used to obtain a solution. 

Using the model, a PAFC system, that was developed under NASA Grant 
NCC 3-17, was improved through the optimization of the heat exchanger 
network. 


Three types of 
evaluated to obtain 


cooling configurations for cell plates were 
the best current density and temperature 
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CHAPTER 1 
INTRODUCTION 

A primary objective of the energy conservation efforts is to reduce 
the consumption of premium fuels such as oil and natural gas. This can 
be realized by displacing the premium fuels with noncritical fuels such 
as coal, or by developing energy conversion systems that can 
siginificantly increase the utilization efficiency of fossil fuels. 
Fuel cells can achieve both objectives by efficiently using a variety of 
fossil fuels, from which hydrogen can be derived, in applications where 
electric power and process heat are simultaneously in demand. In 
addition to their efficiency and fuel flexibility, fuel cells are 
potentially attractive because of modularity, environmental 
acceptability, and adaptability to both utility and on-site integrated 
energy systems. Thus, fuel cell system analysis efforts which are 
essential for design optimization to improve system performance and cost 
effectiveness have acquired significant practical importance. 

For a 'device that was. invented in 1839 by Sir Uilliam Grove, the 
fuel cell has taken a long time to come into its own. It did furnish 
power for the Gemini and Apollo spacecraft and the space shuttle, but 
that was an exotic axid expensive application. Now the fuel cell, vastly 
improved over the spacecraft versions and in assemblies 1,000 times 
larger than the ones carried by the spacecraft, appears to have reached 


1 




, s^age -hare la can «aha a significana conarlbualcn ao a naalcr. s 
supply of electricity. 



3 


1.1 Definition of Fuel Cells 

”A fuel cell is an electrochemical cell which can continuously 
change the chemical energy of a fuel and oxidant to electrical energy by 
a process involving an essentially invariant electrode-electrolyte 
system" (Ref. 1). 

A fuel cell consists of two electrodes, a positive electrode (the 
cathode) and a negative one (the anode), separated by an electrolyte, 
which transmits ions but not electrons. A fuel, typically hydrogen, is 
supplied to the anode and oxygen (in air) is supplied to the cathode. A 
catalyst on the porous anode causes hydrogen molecules (K 2 ) to 
dissociate into hydrogen ions (H+) and electrons. The hydrogen ions 
migrate through electrolyte to the cathode, where they react with 
electrons (supplied through the external-circuit load) and oxygen to 
form water (K 2 O). This electrochemical transformation is isothermal; 
that is, the fuel cell directly uses the available free energy in the 
fuel at its operating temperature. Thus, it is not Carnot-cycle limited 
and can yield a high fuel to DC power conversion efficiency. Unlike a 
battery, a fuel cell does not run down or require recharging, it will 
operate as long as both fuel and oxidant are supplied to the electrodes. 
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1,2 Fuel Cell Efficiency 

The main attraction of fuel cells is efficiency: A fuel cell can 
convert chenical energy directly into electricity and heat. Therefore, 
it is important to mention the theoretical efficiency of fuel cells 
here. The following includes a comparison of maximum fuel cell 
efficiency with maximum theoretical efficiency of conventional power 
generation techniques to illustrate the basic attractiveness of fuel 
cells. 


The process of galvanic oxidation that produces electricity in a 
fuel cell resembles that of a conventional storage battery. However, 
the fuel cell does not store electrical energy; instead, it combines 
electrochemical oxidation (a property of the storage battery) with 
continuous feeding of sufficient fuel (a feature of the mechanical 
engine) to maintain a desired output. 

The efficiency, £, of steam or internal combustion engines, which 
are two of the most widely used conventional methods of generating 
electricity, is limited by the temperatures at which heat is supplied 
(T2C) and rejected (TIC), according to the Carnot cycle- 

c - T2C-T1C 

^ ■ T2C (1-1) 
The maximum practical theoretical efficiency for heat engines is 40% to 
50%. Fuel cells convert energy isothermally (not Carnot-cycle limited), 
and most of the chemical energy of the fuel may be converted to 
electricity. 
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In a H 2 /O 2 fuel cell, the following reactions occur: 

Anode Reaction K 2 = 2H'*’ + 2e" (1-2) 

Cathode Reaction 2e“ + 1/2 O 2 + 2K* = H20(g) (1-3) 

Overall Reaction H 2 + 1/2 O 2 = H20(g) (1-4) 

The overall reaction occured at 177 *C, based on the lower heating 

values (LHV), 

G = “52.9 kcal/g-mol 
H = “58.2 kcal/g“mol 

The maximum amount of heat energy that can be produced isothermally and 
isostatically by the overall reaction is the enthalpy change, A K, but 
the fuel cell can only convert to electricity an amount equivalent to 
the Gibbs free energy, A G. The difference, is the 

oiinimum amount of heat that is produced in fully reversible processes 
both on cathode and anode sides. Therefore, the maximum efficiency of a 
H2/02 fuel cell at 177 "C is 


e 


max 



(100) = 


“52.9 

“58.2 


) 


(100) = 91.0 % 


(1-5) 


The standard potential of a fuel cell, E° , is directly related to 
the Gibbs free-energy change for the overall reaction, according to the 
expression 

Ag = -nFE* - - „ (1-6) 

where n is the number of electrons transferred, or g-equiv/g-mol, and F 
is the Faraxiay constant,, which is equal to. 23,06 kcal/V“g-:equiv. 

In a fuel cell,' all of the f ree 'ener^,'A G, of the reaction is 
available as cell electromotive force (emf) and a large fraction of the 


TAS heat can be recovered in a bottoming cycle or cogeneration scheme* 
In contrast, in any heat engine a smaller fraction of the AH of the 
combustion reaction is available for useful work as a result of the 
Carnot efficiency limitation. 

Nevertheless, there are inefficiencies in a fuel cell, for example, 
the irreversibility which will be discussed in Section 2.2.1, fuel 
utilization, and electronic and ionic conduction resistance through 
components and at contact points. Deviation from ideal behavior caused 
by irreversibility is Illustrated in Figure 1*1. Therefore, under 
realistic operating conditions the efficiency of fuel cell is between 40 
to 45 percent (Ref. 2). 



(Txirrent Density, A/ft of cell 


Figure 1*1 Schematic Fuel Cell Polarization Curve 
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1.3 Types of Fuel Cells 


Fuel cells are usually classified according to the type of 
electrolyte and the operating temperature, as shown in Figure 1-2. h 
brief introduction to the four fuel cell types follows. 


Operating 

Temperature 


Electrolyte 


Electrolyte 

Support 
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(1) Solid Oxide Fuel Cells. Fuel cells with solid oxide electrolytes 
operate at temperatures between 700 *C to 1000 *C using non-noble metal 
electrodes. The electrolyte is usually yttria or calcia stabilized 
zirconia (Zr02) (Ref. 7), which conducts oxygen ions at the operating 
temperature. Solid oxide fuel cells are Intended to offer the advantage 
of burning electrochemically conventional carbonaceous fuels with air, 
without posing any particular electrocatalytical problem. This 
technology is still at the research stage. 

(2) Hoi ten Salt Fuel Cells. This type of cell operates at temperatures 
of 600 *C to 700 *C, using impure hydrogen and air, non-noble atetal 
electrodes, and an electrolyte of molten alkali-metal carbonates in a 
porous ceramic carrier. As with solid oxide cells, high-quality waste 
heat is available for use in cogeneration to increase the overall system 
efficiency. tlolten salt fuel cells are entering the pilot plant 
technology stage. 

(3) Basic Fuel Cells, hqueous KOH electrolyte fuel cells operate at 
lower temperatures (50 to 150 *C) or medium temperatures (210 to 235 
*C>. The latter is called a Bacon cell fdiich was used as the foreruner 
of the Apollo fuel cell system. The alkaline electrolyte used in these 
fuel cells cannot tolerate oxides of carbon in either the fuel or the 
oxidant, therefore, their nonspace use is limited by the cost of 
producing reactants that are free of carbon oxides. 

(4) hcid Fuel Cells. Many acid electrolytes have been considered for 
use in fuel cells because acids are tolerant to carbon dioxide, idiich 
allows the use of lapure hydrogen and air. Phosphoric acid fuel cells 
are the siost advanced of all the fuel cell technologies: Pilot plants up 
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to 4.6-m4 have been successfully operated (Ref. 6) both in New York city 
and Tokyo. Meanwhile, UTC (United Technology Corporation) is developing 
an 11-MI4 plant and l4estinghouse Electric Corporation is working on a 
7.5-MU design, also using phosphoric acid. Because of these facts. This 


I 


research concentrates on this particular system which is called the 
phosphoric acid fuel cell (PAFC) power plant system. 
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;>ower Section Power (Eonditionar Powe 


Spent Fuel 


Figure 1~3 Fuel Cell Power Plant Modules 


Narural gas, methanol, and naphtha are principly considered for 
fuel cell use, so these three fuels are the fuel options treated in the 
present model. Production of hydrogen, thich is the major function of 
fuel processor , occurs by reaction of the fuel with steam. The overall 
reactions are: 


Cn Hm + 2n H 2 O * nC 02 + (2n + 2 ) H 
for naphtha and methane, 'and 


CH 3 OH + H 2 O s CO^ +3% 




for methanol 


(1-7 > 

5 '^^;;; is. 






- • r - 1 - . 


The major components in the fuel processor subsystem are; the 





reformer, two shift converters, and several heat exchangers. The 
reformer Is basically a nonadiabatic, nonisothermal catalytic reactor 
which can operate as high as 1200 *C and 10 atm. The shift converters 
are simulated as adiabatic packed bed catalytic reactors which produce 
hydrogen through a water shift reaction: 

H20 + CO = K2 + C02 (1-9) 

The heat exchangers are modeled as shell and tube type heat exchangers. 

Within the fuel cell stack subsystem (Power Section), hydrogen and 
oxygen react to continuously produce DC electricity, waste heat, and 
steam (as a reaction product). The reacting oxygen is obtained from air 
and the waste heat produced can be removed by passing a coolant through 
the fuel cell stack. For the present simulation, air is used as the 
coolant. 

The reaction which take place at the anode and the cathode, 
respectively, are; 

H2 = 2K^ + 2e- (1-2) 

1/2 02 + 2H* + 2e- = H 2 O (1-3) 

and the overall reaction is 

H 2 + 1/2 ©2 = H 2 O + energy (1-4) 


The electrical power production in the fuel cell stack is 
accompanied by approximately equal amounts of . heat generation which 
should be removed to avoid overheating of the stack. A store detailed 
system flow diagram is shown in Figure 1-4. 



Diagram of Phosphoric Acid Fuel Cell Power Plant System 
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■the catalyst/electrolyte interface, it is electrochemically oxidized 
according to Equation <1-2 ). The electrons are transported through the 
external circuit, and the hydrogen ions are conducted to the cathode 
through the 96% phosphoric acid electrolyte held in place by and inert 
Inorganic or polymeric matrix. 


Balk gas in manifold 

PoroTis conductive backing sheet 
Porous catalyzed layer 

Mateix, H1PO4 electrolyte in wick-type structu 
Porous catalyzed layer 
Current pick-up 


Figure 1-5 Elemental Phosphoric' Acid Fuel Cell Schematic 


The basic cathode structure and materials are the same as the 
anode's (the cathode is not doped) although the physical parameters 
(pore size, thickness, etc.) may be different to achieve the desired 
mass transfer characteristics. Hhen the hydrogen ions arrive at the 
cathode reaction sites, they react with the oxygen, which has diffused 
through the pores of the cathode's Teflon region. The reaction is 
Equation (1-3). 


Single cell assemblies are stacked in a series-connected bipolar 
Bode. The directions of flow channels for air and hydrogen fuel are 
perpendicular to each other. In addition to the DC power generated, 
heat is produced as a result of the reaction. This heat is removed from 
the cell stack by cooling air passing through channels in cooling plates 

OF 
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located approximately every fifth cell so as to maintain the cells at 
the desired temperature. Figure 2-17 is one of the symmetric parts in 
the stack. The heat carried away by the cooling air is used to generate 
low pressure steam to drive a steam turbine. The cooling air is then 
recirculated back to the fuel cell system. 

The effluent fuel from anode side goes to the burner after mixing 
with air from the air compressor. The cathode exhaust is first allowed 
to pass through an expander to recover some power. It is than cooled 
to condense out water in the stream. The condensed water is separated 
from the gas stream and sent to storage so that it can be used for the 
fuel processor or steam turbine. 

1.41.3 Power Conditioning Subsystem 

The primary functions of power conditioning subsystem are to 
control the magintude of the direct current DC generated in the fuel 
cell stack and to combine and convert the DC to alternating current 
(AC). Secondary functions are to povide appropriate fault protection, 
reduce harmonic distoration, and to provide power factor correction 
capadsility. The DC power source will be connected in various 
series/parallel arrangements to deliver the power at the required 
voltage and current. 


CHAPTER 2 


DESIGN OF PAFC SYSTEM 

A primary objective of the energy conservation efforts is to reduce 
the consumption of premium fuels such as oil and natural gas. This can 
be realized by displacing the premium fuels with less critical fuels 
such as coal* or by developing energy conversion systems that can 
significantly increase the utilization efficiency of fossil fuels. Fuel 
cells can achieve both alternatives by more efficiently using a variety 
of fossil fuels* from which hydrogen can be derived, in applications 
where electric poi/er and process heat are simultaneously in demand. 
However* while the potential conservation of fuel is important* the 
ability of fuel cell to penetrate the marketplace is essential in order 
to succesfully achieve that potential. Thus* fuel cell system analysis 
efforts* which are essential for design optimization to improve system 
performance and cost effectiveness have acquired significant practical 
importance. 

- ^ The purpose of this chapter is to utilize mass and energy balance 
techniques with electrochemical and chemical Idnetic disciplines to 
develop mathematical models and the associated digital .! .computer; 
simulations. The purpose of tlie coiqputer simulation is to develop. a 
simple but sufficiently accurate model of a PAFC syst«n to allow for 
careful engineering evaluation of the system. The temperature* 
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pressure, and naterial balance at various locations in the systen and 
the sensitivity of operating conditions in major components are 
determined using the computer code. In addition, alternative design 
configurations in the reformer and fuel cell stack will also be examined 
for the optimal design of the system. The philosophy of the computer 
simulation is not to develop a highly accurate, rigorous mathematical 
analysis in microscopic detail but rather a workable but sufficiently 
accurate approach to a practical engineering problem. 


In this chapter we consider two kinds of mathematical models for 
the components. One is the lumped model and the other is the 
distributed model. The lumped or "single-staged" model calculates the 
material and energy balances and the physical and chemical properties of 
the streams at one operating condition, such as the mean temperature and 
pressure. On the contrary, the distributed or "multi -staged" model, 
which usually is solved by the finite-differences method, handles the 
balances and properties at different sections. In the reformer and both 
shift converters, the lumped model utilizes the approach differential 
temperature (ADT)*» however, chemical kinetic will be considered in the 
distributed model . 


■ There are two schemes of PAFC system idiich will be simulated in 
this chaprter for steady state performance. One is the design develop)^ 


*: The approach differential temperature is the temperature difference 
between the exit temperature and the assumed operating temperature at 
t^ch equilibrium occurs. 
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with this research which will be referred to as the C$U design, the 
other is the Uestlnghouse conceptual design (Ref. 12) which was 
described in Chapter 1. The CSU design scheme can be described as 
follows: 

As shown in Figure 2-1, methane which is circulated by compressor 
<C) is preheated by heat exchanger E-1 prior to mixing it with the super 
heated steam which receives its heat by passing through heat exchanger 
E-9. Before entering the reformer, the methane steam mixture is heated 
via heat exchangers E-2 and E-3. Inside the reformer, methane is 
catalytically reformed by reaction with excess steam to produce carbon 
monoxide, carbon dioxide, and the desired product, hydrogen. The 
effluent from the reformer is cooled by flowing through heat exchanger 
E-2 before it enters the high temperature shift converter S-1. The 
function of the high temperature shift converter is to increase the 
hydrogen concentration and to reduce the carbon monoxide concentration 
of the reformer gas effluent. The temperature of the effluent from the 
shift converter S-1 is then reduced by passing through heat exchangers 
E-1, E-9, and E-6 before entering the lew temperature shift converter 
S-2. The low temperature shift converter further increases the hydrogen 
concentration by promoting the shift reaction at a lower operating 
temperature. The effluent from the low temperature shift converter then 
enters the fuel cell containing H 2 , CO, CH^, CO 2 and K 2 O. The fuel cell 
converts inputs of hydrogen and oxygen to DO power, water and . heat. 
Oxygen is delivered to the fuel cell by air compressor A which " also 
provides air to the reformer burner. The spent fuel from the fuel cell 
anode goes to the burner after mixing with air supplied' by compressor A. 




















Before entering the burner, the Biixture is preheated by the burner 
effluent via heat exchanger E-<1. The spent fuel is then burned vith 
whatever additional methane is needed to provide the thermal energy* 
necessary for the reformer reaction. 
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Heat generated in the fuel cell is removed by heat exchangers E-7 
and E-10. Heat from heat exchanger E-7 can then be utilized in 
Industrial heat processing or space heating and cooling, while exchanger 
E-10 is used to preheat the water supplied by liquid separator Q to 
provide the necessary steam needed for the reforming process. The 
effluents from the burner and fuel cell cathode will have their water 
removed and separated by condenser E-5 and liquid separator Q before 
allowing them to be exhausted to the atmosphere. 

This chapter will begin by modeling the components in the fuel 
processing subsystem and fuel cell subsystem, then the system steady 
state performance will be presented. 
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2.1 Hodeling of Fuel Processing Subsygten 


Production of hydrogen, which Is the major function of the fuel 
processing subsystem, occurs fay reaction of the fuel with steam. The 


temperature shift converter, the low temperature shift converter, and 
several heat exchangers. 


2.1.1 Heat Exchaneer 


A zero capacitance sensible heat exchanger is modeled in the 
double-pipe counter mode. 


Lumoed hodel 


For the counter siode, given the hot and cold side inlet 

temperatures and flow rates, the effectiveness is calculated fcr a given 

fixed value of the overall heat transfer coefficient. The mathematical 
■ 

description which follows is covered in detail in Ref. 20. 


Tbo - Ihi - E ( ) (Thi - Tci) 


- E (-^^) (IM - Iti) + Tel 


» ECmin (Thi - Tci) 


C m ■in 


[ 1 — Cmija/Cnax ) 


1~CC m<n/ C mtt-r )e ^'®in ■ : 
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where Cc; capaci^y rate of fluid on cold side, BcCpc, J/s-K 
Ch: capacity rate of fluid on hot side, MhCpc, J/s-K 
Cmax: naximum capacity rate, J/s-K 
Cmin: ninimum capacity rate, J/s-K 
Cpc: specific heat of cold side fluid, J/g-K 
Cph: specific heat of hot side fluid, J/g-K 
E: heat exchanger effectiveness 
lie: fluid mass flow rate on cold side, g/s 
ilh: fluid mass flow rate on hot side, g/s 
Qt! total heat transfer rate across heat exchanger, J/s 
Tci: cold side inlet temperature, K 
Tco: cold side outlet temperature, K 
Thi: hot side inlet temperature, K 
Tho: hot side outlet temperature, K 

DA: overall heat transfer coefficient of exchanger, J/m*-s-K 
Distributed llodel 


From the first law - of thermodynamics, the energy balances for 

tubeside gas, shellside gas, and tubewall can be written as, 

T, ,, .-TTurr ,1 Pi n (d:T Qi dt 

Tubeside: hi 7t di CTw-t) = j. 

T 

Shellside: ho do <T-Tw) = -^o fto Cpo Vo - 
Tubewall: ho do (T-Tw) - hi di (Tw-t) =0 
where h: composite heat transfer coefficient, J/sec-m^-K 

density, g-mole/m* ,, ^ . . u., . > 

Cp: heat capacity, J/g-mole-K .^ 
t: temperature of tubeside gas, K 


t' 


H 


T: temperature of shellside gas, K 


Tw: temperature of tubewall, K 


V: velocity, a/ sec 


d: diameter, m 


subscript 1: tubeside 


o: shellside 


w: tubewall 


Boundary conditions 


zsO t=t inlet 


at=L T=T inlet 


These simultaneous ordinary differential equations were solved by 


finite~differehce method. The diagram of one section is shown in Figure 


VoJ.Ti 


■Voj+2,Tj+2 shell side 


tubeside 







Figure 2-2 Secticn J of heat exchanger 


Itesults 


Distributed model was used to estimate the dimensions of the heat 


in Figure 1-3, Table 2-1 shows the results where the known 


inlet and outlet temperatures are from the Uestinghouse conceptual 
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.} 

design (Ref. 12 ). These results will be used as input data for the 
system steady state and transient state simulations. 

2.1.2 Shift Converters 

The function of both types of shift converters (high temperature 
and low temperature) Is to further increase the hydrogen concentration 
and to reduce the carbon monoxide concentration of the reformer gas 
effluent. The equation, 

CO + HgO = H 2 + CO 2 (water shift reaction), 

dominates the Biaterial changes in the shift converters. The methanol 
input fuel does not need to pass through shift converters because the 
carbon monoxide level is low. 

Lumoed hod el 


In this model, the water shift reaction is assumed to be at 
equilibrium and utlizes an approach differential temperature (RDT) 
ranging from 15 to 25 *F. The material balance is 


- PcOj Phj (FcOj-*OC) (FHi 4-x) 

fco <Fco-x) (Fh^o-x) (2-1-1) 

where Kz: equilibrium constant of shift reaction at ftDT 
^ P : partial pressure of component, atm 

F : inlet molar flow rate of component, g-mole/s 
X : reacted amount rate, g-mol^s 

Equation 2-1-1 can be solved for x. Newton's method was used in the 
computer program.. 
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The energy balance equation for the gases in the shift converter 


includes the reaction and sensible enthalpies. For adiabatic the 


process in the shift converter 

nj <Ah“f)j - ZZ Txi (^h‘m + Z nj 

P5 R5 P5 

- 21 ni ) <Cp)i dT = 0 (2-1-2) 

where the subscripts PS, RS correspond to the products and reactants in 


the shift converter, respectively. Tf, Ti are the final and initial 


temperatures of the gases, respectively. The only unknown in the 


equation ,Tf, is determined iteratively. 


r 

(Cp)j dT 

J29& 


The Ergun equation, which estimates pressure drop caused by the 
flow of gas through dry packings, is used to determine the pressure drop 
in shift converter and reformer. The equation is (Ref. 22) 


1-6 G / " 150 g-€) M l \ . 
Ap = 1878* C^-dpi^pl^ dpG +1.75/h 

where £ : void fraction in bed 

J { : viscosity, Kg/m-s 

dp: effective diameter of packing particle, m 
G : superficial gas mass velocity, Kg/s-m 
h : packed height, m 
P : density, Xg/m® 

pressure drop, atm . 

Distributed Model 


(2-1-3) 


In the distributed model simulation, . the „^^thift^ convert^ , is , 
simulated as an adiabatic packed bed catalytic reactor tdilch produces 
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hydrogen through the water shift reaction. 

The general assumptions made for simulation of a packed bed reactor 
(like the reformer and shift converter) are discussed here; 

(1) Axial dispersion and radial gradient are negligible - plug flow 
condition. Generally if the ratio of the length of the 
reactor to the catalyst's diameter is greater than 100, 

the axial dispersion effect is negligible. 

(2) A uniform temperature exists throughout each catalyst particle, 
and this temperature is the same as the gas temperature in that 
section of catalyst bed. 

(3) The kinetic expression represents a global rate, and therefore 
neglects reactivity differences found between the inside 

and outside of the catalyst particles. 

(4) Entrance effects are negligible. 

(5) Heat transfer by radiation is negligible. 

(6) A single reactor tube is analyzed. Thus, all the tubes in the 
reformer or shift converters behave independently of one another. 

(7) Ideal gas behavior is assumed. 

(8) The outside shell wall is adiabatic. 

A more detailed discussion of assumptions (2) and (3) is provided in 
~ Appendix 2 by examination jof the "internal" and "external" effectiveness 
factors of commercial catalysts used in reformer and shift converters. 

Mass Balance 



From the general continuity and assun^ions, the kinetic mass 
balance for the shift converter is 
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d2~ e 

where V : average velocity of fluid through the bed, si/s 
C : g-mole of CO per fluid 
ra' : reaction rate, g-mole of CO/s-Kg catalyst 
e^: density of catalyst, )Cg/m^ bed 

Reaction Rate 


(2-1-4) 


From Ref. 24 the rate expression for a high temperature shift 
converter with a catalyst composed of Fe20^ and Cr20^ is (Ref. 23) 

V- - It^PHaO Pco -Ka Pcoa PHa) 

' ~ A PhjO + Pco;i (2-1-5) 

where r : reaction rate, m*C0 3 STP/m’ bed-s 
k : rate constant 
A : constant 

for a low temperature shift converter with a catalyst composed of ZnO, 

Cr 203 , and CuO, the first order rate and Arrhenius expression is 

-EA/flT 

-r = Ko e Pco (2-1-6) 

where r in g-mole of CO/s-kg cat. 

Ko: Arrhenius frequency factor, g-mol/s-kg cat. -atm 
EA: activation energy, J/g-mol 
R : gas constant 
T : temperature, K 

Energy - Balance ^ ^ 

si I - The energy .. balance ^f or the reacting gas is 

Cp p 


(2-1-7) 
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„Here Cp: hea-- capacity =5 fluid, .i/g-cola-n 
p: fluid density, g-mol/m^ 

-AH: heat of reaction, J/g-mol of CO 
^ : temperature of reacting gas 


These two oridinary 
finite-difference method. 


differential equations were solvec ty 
A diagram of one section is shown in 


the 

gure 


2 - 3 . 



Plgur. 2-5 S«tlcn 3 of shift converter 


Results 

The distributed model was used to estimate the dimension- 
high temperature anh_ low temperature shift converters in Figure 1-3. 

w»wn in Table 2-2, where the known inlet and outlet 
The results are shown in laoxe 

_ 12 These results will be 

caaposltlons and temperatures are from Ref. 12. Th 

a « for the steady state a«i transient response 
used as input data for tne sxea y 


simulations of the system. 


:-agh Temperature shift Converter 
Islet: 


\now Rate 
Flow So^'\^ 

CH4 

CO 

C02 

820 H2 N2 

02 

m 

1 ' ! 

P 

1 

AU 

11.78 

100.16 

44.6 

205.66 472.96 7.57 

0 

I 720 

1 5.1 


Units: 51ov Hate: Ih-aole/hr 

Temperature: ? 

Pressure: atm 


Outlet: 

■ ,?lou Hate 


Hesouroe- 

0H4 

CO 

C02 

H20 

H2 

K2 

02 

T 

P 

Westinghouse 

Simulation 

11.78 

48.2 

96.56 

151.7 

524.92 

7.57 

0 

840 

1 

4.8 

Our Results* 

11.76 

49.8 

95.0 

155.5 

525.4 

7.57 

0 

864 

5.0 


• A single Pipe reactor is oonsideree with these conceptual design parameters: 


SC. of reactor tubes 92 

ID of reactor tube 0.5 ft 

Height of reactor tube 6 ft 


I<ov Temperature Shift Converter 
Inlet: 


Flow Hate 


F1 ow' 1 jo>\^ 

CH4 

CO 

C02 

H20 

H2 

K2 

02 

T 

P 

ai 6 

! 11.78 

i 

46.2 

96.56 

151.7 

524.92 

7.57 

0 

597 

4.42 

Cutlet; 

>^low Hate 1 
3esour5>^s^^ 

CB4 

CO 

C02 

H20 

K2 

S2 

02 

I T 

? 

.estin^ouse 

Simulation 

11.78 

10.8 

155.96 

114.5 

562.5 

7.57 

0 

495 

i.42 

Our Hesulte* 

11.78 

10.9 

153.8 

114.4 

562.2 

7.57 

0 

502.2 

4.56 


* NO. of reactor tubes 
to of reactor tube 
Height of reactor tube 


92 

0.25 ft 
2 ft 


Table 2-2 ■Estimation of dimensions of shift converte 
in Westinghouse FAFC system 


OF POOHijUALITY 
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2.1.3 Rgfonrier 


Tbe key component in the fuel processing subsystem is the reformer 
which catalytically reforms methane (methanol or naphtha) by reaction 
with excess steam to produce carbon monoxide, carbon dioxide, and the 
desired product, hydrogen. Methane will be the only input fuel in the 
following discussion. 

Two reactions are assumed to be the principle reforming reactions 
in the model, they are 

CH 4 + H;jO = CO + 3Hj (demethanation reaction) 

CO + HjO s CO 2 . + Hg (water shift reaction). 

Appendix 1 lists all of the possible reactions and discusses the minimum 
steam to carbon ratio (S/C) required to avoid carbon formation. 

Lumped Model 

In the lumped oiodel both of the reactions, demethanation and shift 
reaction, were assumed to be at equilibrium by utilizing the respective 
ADT's of each. The equilibrium constants were determined from the 
temperature. The equilibrium expression are 

K - PcOa _ Ycoa Y^a 

•^1 *" PCH 4 pHaO Ychu, YH 20 (demethanation) 

1/ ~ Pcoa Phz _ YcoxYHz 

•^2 Pco PHzP Yco Yh^o - . (water shift), 

vhere K 1 anid K 2 are the equilibrium constants of demethanation and water 
shift reaction, respectively. Expressing the mole - fractions as the 
individual molar flows divided by the total molar flows yields: 


3 ^ 


and 


u ( Fco-x^y) (Fh 2 +x+3y)^ 

(FcM^-y) (FWaO-x-y) (Fr+2y>* 

1/ - (Fco;+x) (FH2-*^~*'3y) 

^2- (Fco-x+y) (Fh^ o-x-y) 


( 2 - 1 - 8 ) 


(2-1-9) 

where y is the conversion amount rate in the demethanation reaction and 
F is the total inlet flow rate. Equations (2-1-8) and (2-1-9) can be 
solved for x and y. Newton's method was used in the computer prosran. 
Figure 2-4 and 2-5 show the equilibrium conversion of methane at 
different S/C ratios, temperatures, and pressures. It can be observed 
that greater pressures and greater S/C ratios will favor the conversion, 
but over the temperature ranges, temperature seems to be the major 
factor. 


The quantities involved in the energy balance will be the sensible 
enthalpies of the gases, the reaction enthalpies of the gases, and the 
heat transferred from the combustion gases to the reformer gases, 

The value of Qg./q can be determined from 


Q 0 .pj= UAAT 7 ft= Hout - Hin 
where, A T^ is the log mean temperature defined as 
(T^c - ? - (T^- 


( 2 - 1 - 10 ) 


^ — 


Jl'n- 




( 2 - 1 - 11 ) 

where, T^ is the temperature of the combustion gases after leaving the 
reformer; T^^ and Tf^^are the temperatures of the reformer gases before 
entering and after leaving the reformer; h is the heat transfer area; 
and U is a modified form of a heat transfer coefficient. 







Thus, from the first law of thermodynamics and equation (2-1-11)/ 


the energy balance for the reformer gases can be written as. 




= XI m. (Ah?)- - ZIn. (Ah?)- 

PR ^ ' <r PR ^ ‘ pr^ 


rliR 


2S6 


(Cpy dT 


- E 




(Cp)j dT, ... 


( 2 - 1 - 12 ) 


rR -^29d 

where the subcripts PR and rR stand for products and reactants in the 
reformer, respectively. 


Distributed Hodel 


Kinetical analysis was used for simulation of the performance of 
the reformer. The reformer is basically a nonadiabatic, nonisothermal 
catalytic reactor, that is heated on the shell side by combustion gases 
from burner. Figure 2-6 shows its simplified scheme. 


Refarner 



Combustion 





di z 


Reforming 


_Pigure_^-6 


Simplified reformer diagram 
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The demethanation reaction is assumed to be kinetlcally controlled 
and, hence, occurs at a finite rate, while the water gas shift reaction 
is assumed to be equilibrium controlled. The modeling assumptions 
discussed in Section 2.1.2 are also applied here. The demethanation 
reaction used in this model is slightly modified with linear 
combinations of the original demethanation reaction and shift reaction, 
which results in 

CH 4 , + 2 H 2 O = CO 2 + 4H2 (2-1-13) 
In the equilibrium calculations, the demethanation reaction choice 
causes no changes in the final results. However, the kinetic 
consideration will cause the final results to vary slightly with the 
reaction choice. 

Hass Balance 


From the generalized continuity and the assumptions, the kinetic 


mass balance is 

ra€8 

cLZ ~ ~ £ ' 

this is the same as Equ. 


( 2 - 1 - 1 ^) 

2-1-4 except that the basis is changed to 


methane. 


Various kinetic expressions for the reforming of methane with steam 

have been proposed which could provide the rate equation (Ref. 24, 25, 

2S), The simplest form among the proposed expressions is the first 

order rate expression, which is 

: -eA/«T . , ^ / 

-ra' s Ko e PcH<f (2-1-15) 

in Arrhenius form (the same definitions as Equ. 2-1-6). 


Unfortunately, 
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little agreement can be found for the values of the kinetic parameters, 
some values may be three orders of magnitude different from others. The 
data from Ref. 26 using a commercial catalyst (Gindler G-5SE) is used in 
this model. 

The vater gas shift reaction is assumed to be at equilibrium. The 
conversion quantity is based upon the carbon dioxide mass balance. 
Thus, when coupled with the demethanation reaction, the water gas shift 
reaction proceeds in reverses, therefore, the shift conversion is always 
negative. Using these two reaction schemes, all of the molar flows 
anywhere in the reformer can be written in terms of the feed quantities 
and the conversions of the two reactions. 

Energy Balance 

Two energy balances are required for the system: one for the 
reformer gases and one for the combustion gases. The reformer gas 
balance includes its own sensible heat change, reaction enthalpies, and 
heat transfer from the hotter combustion gases. The combustion gas 
balance involves sensible heat change and heat transfer. This 
translates quantitatively into Equation (2-1-16) and (2-1-17). 

pMV Cp-^= <-AH,)-^+ hi;tdi (Tw-t) (2-1-16) 

9o Vo Ao Cpo-^= -hoXdo (T-Tw) (2-1-17) 

where AH,: demethanation reaction enthalpy, J/g-mol CHt,. ■ 

AH^: water shift reaction enthalpy, J/g-mol CO ■ 

M : inner tube cross area, 

hi : wall heat transfer coefficient of tube side, J/s-m*-K 
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Tw : wall reinpera'ture,K 

di : inner tube diameter, m 

T : combustion gas temperature, K 


do ; outer tube diameter, m 


subscript o refers to the combustion gas side 


1 


m 


There is greater uncertainty in estimating the heat transfer 
coefficient at the wall of tube than the rate expression. The scatter 
in experimental data is very high (Ref. 22, 24, 25). The situation will 
be even more complicated by considering the unequal stoichimetric 
reaction (Ref. 28). Due to Beek's recommendation (Ref .29) the modified 
Thoenes-Kramers (Ref. 30) correlation should be used for sphere-like 
particles near the wall, which are used in the model. 


hi(dp/kjT) = 2.58(Re)^ (Pr)^+ .094 (Re)®* (Pr>®** 
where dp: equivalent particle diameter, m 
k^: thermal conductivity, J/s-m-K 
Pr: Prandtl number 


(2-1-18) 


Re: partical Reynolds number 


m 


F- 




Differential equations, (2-1-1), (2-1-14), (2-1-16), and (2-1-17) 
were solved simultaneously with the inlet conditions as the boundary 
conditions. The Ergun equation (Equation 2-1-3) is used to evaluate the 
pressure drop. The Computer code which was developed to simulate the 


performance used the finite-difference method to solve these 




n 



simultaneous ordinary differential equations. Figure 2-7 shows one 
section of the reformer, and the program flowchart for the model is 


given in Figure 2-8. 
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Results and Comparison 

The model results were compared to the Phillips Petroleum reformer 
data (Ref. 25). The reformer described in Ref. 25 is used in the 
refinery industry, so the tube is electrically heated to maintain a 
vertical temperature profile. The model discussed in this section were 
modified to fit this special condition. Figure 2-9 gives the 
temperature, pressure, and conversion profiles along the reformer tube. 
Table 2-3 shows the comparison together with Mestixighouse BOLTAR 
simulation results (Ref. 11). As the results indicate, the model 
corresponds to the given data very well. Pressure drop considerations, 
different rate expressions, and different heat transfer coefficient 
estimation enable the CSU model to better simulate the Phillips 
Petroleum reformer in comparison to BOLTAR. 

The CSU model can also be modified to simulate the "regenerative" 
reformer which will be discussed in Section 2.2.3. Because the 
reformer's performance greatly depends on the temperature level, the 
unceiT;ainty in the rate expression and heat transfer coefficient will 
limit the simulation model accuracy. 




Position Batio Along Refaraing Gas Plow 
^gure 2-9 The reforming tube profile obtained for the distributed model 





Reformer Results Comparison 
Phillips Petrolexim 


Input Data 


Reformer Data 
(Ref. 25). 


BOLTAR(Ref. 13 ) CSD 


Height, m 

12.2 

12.2 

12.2 

Tube ID, cm 

13 . 

D2=13. 

D2=13. 

Tube wall thiclcne3s,cm 

1.35 

1.3 

1.35 

Catalyst part, size, cm 

1.59x1.19x.48 

1.27 

1.27 


Tings 

spheres 

spheres 

Mass flow rate 

26700 

P0=19.1 

P0=19.1 

(superficial) 

kg/hr-m'^ 

kg-mole/hr 

kg-mole/hr 

Steam/methane ratio 

7:1 

7:1 

7:1 

Gas inlet temp., “C 

364. 

364 . 

364 . 

Tube wall temp., «C 

704+31xz-1.2xz2 

same 

same 


(z in m) 



Inlet pressure, MPa a. 

14.3 

13.2 throughout 

14.3 

Output Results 

Given 

BOLTAR 

CSU 

Outlet pressure, KPa a. 

12.2 

15.2 

12.7 

Outlet temp., *C 

793 . 

CM 

CO 

776 

Conversion, % 

91.7 

06.1 

91.9 



Table 2 - 3 Con^arison of CSU model with experimental 

data and BOLTAR 



OF POdK 
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2.2 Modeling of Fuel Cell Stack Subsystea 


In the fuel cell power section, air, in excess of the 


effluents from the low temperature shift converter enter at the anode. 
The anode input contains CH^, H 2 O, H 2 , CO and CO 2 . In this analysis, it 
is assumed that a fixed percentage of hydrogen is consumed at the anode, 
and the K 0 being formed exits the fuel cell, with the depleted air, 
through the cathode exit. The overall reaction in the fuel cell power 


section is. 


+ 1/2 Oi = H 2 O 


( 2 - 2 - 1 ) 


This section presents the two distinct mathematical models of fuel 
cells that were developed and computer programs that were written to 
perform the necessary calculations. 


The simplified lumped model is an "input-output” model developed 
for the system trade-off studies. 


The detailed distributed model is a finite -difference model of the 
operation of the fuel cell which was used to calculate the effects of 
the cell and module design on performance. It calculates the current 
density distribution in the cells as a function of the local reactant 
compositions, focal temperatures, 'catalyst utilization factors, etc. 
Since these are interdependent (e.g. the .local traperature depends on 

‘ ^ ‘ f'* f 

the local current’ density), the computations are highly iterative and 
require considerably more computer capacity and tine than the lumped 
Bodel. An associated computer program will be used to compare an 


alternative design of cooling scheme in the stack. 


2.2.1 Lumped Model and Voltage-Current Characteristics 


The simplified lumped model provides a rapid (in terns of computer 
time) means of calculating the fuel cell nodule output characteristics 
(voltage, current, and heat generation rate) in terms of the inputs from 
the fuel processing subsystem and the gross fuel cell design parameters 
such as catalyst loading. 


The mass balances of hydrogen, oxygen and water are as follows; 






KXh 2 ~ - (Imean • A)/(n ?> 

NXq^ = NIoi - (Imean • A)/(2n 
KXhjO “ NIh^o (Imean • A)/( n yo 


( 2 - 2 - 2 ) 


(2-2-3) 


(2-2-41) 


where N){; exit flow rate of hydrogen, oxygen, or steam, g-mole/sec 
NI; inlet flow rate of hydrogen, oxygen, or steam, g-mole/sec 
Imean: mean current density, h/cm^ 

A: effective area of cell plate, cm^ 
n: number of Faraday equivalents transfered 
Faraday constant 


The energy balance for the fuel cell is, 

-<8 . U.) i ZZ (A h°f ): - H h f 

. .PP ^ f hF ^ 


(2-2-5) 





'*ere the subscripts PF, irF represent the products and reactants' in the 
fuel cell, respectively. Tfp is the final temperature of the products 
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and TiF is the initial tenperature of the reactants in the fuel ceil. 
The nj and n.i are the species flow rates of the products and reactants, 
respectively. The terns Q and H are the rates of heat and the 
electrical energy generation by the fuel cell, respectively. Q is 
proportional to the specific heat generation Qp where: 

Q = NpXnTnQp (2-2-6) 

and Qf = (^*' - V) I (2-2-7) 

where Q : total heat generated, J/sec. 

Qp : heat generated per unit area of cell, J/sec. cn* 

Np : number of cells 
Xn : width of cell plate 
Yn : length of cell plate 
I : fuel cell current density. A/cm* 
s heat of reaction, J/g-mole of 

Voltage-Current Characteristics 

Because of the irreversibility, the voltage V for a working fuel 
cell is the difference between the open circuit voltage ,and the cell 
polari 2 ation terms: 

V - Z (2-2-8) 

*^h®re E: Nernst potential (reversible open circurt E.H.F.) 
overpotential or jpolarization 

The reversible cell potential, E is given by xHe Nerh« equation-'-- 

c- JRT f YHaV PtYOa. \ 

£<,= E (T) + ZF In V YHaO ' ' ' ' '''' ^ ‘^'- ^ (2-2-9) 

Pt: total pressure, atm 


qS 


Eg (7): standard E.K.F. of cell at temperature T, volts 
Eo<T)= 1.261-0.00025 7, T, K (Ref. 10 ) 

YH;j : mean mole fraction of hydrogen at anode 
: mean mole fraction of oxygen at cathode 
YH 2 O: mean mole fraction of water vapor at cathode 


The polarization term "5^ consists of four components, 

7^ = 7(_a +^r +^d +7^co (2-2-10) 

where 7[a: activiation poparization at cathode, volts 
T^r: resistance polarization, volts 
7|^d: diffusion polarization, volts 

t^o: activiation polarization at anode due to co poisoning 
of catalyst, volts 


RT i 

=o4oZF In (io)(SR)(CLKCU) 

with c^o: transfer coefficient 

i : current density, nA/cra* 

io: exchange current density of cathode, mh/cm^ 
SA: specific catalyst surface area, cm* /g 
CL: catalyst loading on cathode, g/cm* 

CU: catalyst utilization factor. 


( 2 - 2 - 11 ) 




The exchange current is a function of the acid concentration/ 
temperature, and partial pressure of the oxygen. The 'acid concentration 
Is a function of the waiter vapor partial' pressure' which permits 
^rrelation of io as a function of 'Y02, YH2d, and T. 'empirical fit 


j 
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io = 222.7 (PtYC2 )«»•• . (PtYK20)®**’’’ exp (-6S32/T) (2-2-12) 

The resisTance polarization is 

7[r = ir 

where r: specific cell resistance, ohm-cn* . 

The expression of 7^ co was chosen to have strona; temperature 
dependence, be directly proportional to Yco, and have a logarithmic 
dependence on i, iao and catalyst effective area. The resulting 
expression (Ref. 10 ) is 

co=0.0782Pt Yco exp 9190^“ In (cLa SA CU (2-2-13) 

where CLa: anode catalyst loading, mg 

iao: anode exchange current, mA/cm* 

Diffusion polarization has been neglected here because it is 
significant only at very high current densities. 

In the associated computer program. Subroutine VI, calculates cell 
voltage as a function of the current density or alterTiatively solves the 
nonlinear equation to evaluate current density as a function of the cell 
voltage. 

It is difficult to compare the predictions of a cell voltage 
«<luation with experimental data because of the generally incomplete 
•l®®trode or other cell specifications. Figure 2-10 illustrates the 
cf electrode performance on cell potential. Each of the six 
curves was developed for different catalyst utilization factors, using 




Current Density 
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Equation 2-2-8. The curves are compared with the prediction of the 
design relations from Ref. 10 . The latter are in the nature of 
off-design corrections to a baseline fuel cell operating on pure 
hydrogen and oxygen, 1 atm, pressure, 175*C, and low utilization. 

2.2.2 Current Density Distribution 

In the fuel cell module, the combined modeling of temperature and 
current distribution is an absolute condition for reliable scaling-up of 
the results obtained with small cells, and for predictive models 
starting from elementary porous-electrode representations. 

This subsection describes the calculation of the current density 
distribution over a cell plate on which the air and fuel flows are at 
right angles. The procedure divides a cell plate into "grids" which are 
small enough so that variations in fuel and oxidant composition cund 
temperature are negligible. Then by means of calculation of the 
boundary conditions for each "grid" and iteration, a solution will be 
obtained that satisfies the input specifications (e.g. average current 
density, fuel and air utilization, and reactant flow rates). A diagram 
of the "grid" is shown in Figure 2-11 . 

The overall method is to first specify a desired average current 
density i for the whole plate and then determine the corresponding 
voltage V for the plate. This voltage will be determined such that it 
produces unique local current densities over the plate whose average 
value approximates i within a specified tolerance. A trial -and-error 
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procedure is used to estimate the local current density and overall 
voltage. The model basically applies the same voltage-current equation 
used in the lumped model (described in Section 2.2.1) to each grid 
section of the cell. 

Kathematical Formulation 

Exit flow of hydrogen from grid (i» j) 

= NIh 2<1»J> - <Ki,J) A>/(n F> (2-2-14) 

Exit flow of oxygen from grid <i, j) 

= NIc>i(i,i) - <I<i,J) A)/(2 n F) (2-2-15) 

Exit flow of water from grid (i,j) 

= KI„,o + (I(i,j) R)/( n F) (2-2-16) 

where NXh*. Oi.Hjo * hydrogen (oxygen or water) portion 

flow rate at exit of grid (i,J>, 
g-Bole/sec. 

O 2 H 2 O * hydrogen (oxygen or water) portion 

flow rate at inlet side of grid (i,J) 
g-Bole/sec. 

I (i,j) : current density of grid (i^J), A/cm* 

A : area of grid, cb* 

The flow chart of executive program for calculating current density 
distribution is shown in Figure 2-12 • 

" Results 

A current density distribution over an isbthanul''TOll plate whose 
dimensions are 0.30 m x 0.43 m and vhose average operating tenperature 
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is 450 “K has been detercined from the computer prograr.. This 
distribution is presented in Figure 2-13 as a field plot. of lines of 
constant current density. The fuel and process air conditions are given 
in the figure: where represent utilization of fuel and oxygen 
respectively. The naxinum current density occurs at the corner of the 
plate common to the two inlet streams and is 30% above the average 
current density, however the minimum current density occurs in the 
opposite corner and is 20% below the average. Since there is no 
satisfactory experimental method of measuring the local current density 
(Kef. 11), confirmation of the predicted current density must be 
inferred from temperature measurement. 


The computer program can be modified to simulate different types of 
cell plates - i.e. Zee plate (Ref. 11) and hexagonal type plates. The 
operating conditions and results for these two alternatives are shown in 
Figure 2-14 and Figure 2-15 . 


2.2.3 Thermal Analysis and Temperature Distribution 


The electrical energy production in phosphoric acid fuel cells is 
accompanied by approximately equal amounts of heat energy generation. 
Removal of this heat can be accomplished by a suitable flow of input 
gases or by using separate cooling plates. 


The work reported in this section is directed towards estioiating 
steady state temperautre profiles in practical phosphoric acid fuel 
stacks. The fuel cell stack considered in this section is composed 
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of cell plates on which the air (oxygen) and fuel (hydrogen) flows are 
at right angles. A cooling plate is placed between individual groups of 
cells at a regular interval. Symmetry in the stacking direction occurs 
at the Biddle of a cooling plate and midway between cooling plates. 

2. 2. 3.1 Previous Uork 

Estimation of the temperature profiles in an operating cell is 
important for the estimation of the power density distribution* thermal 
stability* and cooling requirements. Only a limited amount of 
information on this subject has been reported in the past. Baker and 
co-workers recognized this need and have performed a comprehensive study 
of steady state heat transfer in electrochemical systems (Ref. 6*14,15 
). They studied various cases involving one dimensional axialysis of a 
single adiabatic fuel cell and a three dimensional analysis of a 
multicell stack. 

A single fuel cell with no lateral heat transfer and no conduction 
of heat through the cell in the direction perpendicular to the gas flow 
was considered (Ref. 14 ). Heat transfer by conduction in the 
direction of the gas flow was considered negligible in comparison to the 
beat transfer by convection* and axialytical expressions for the 
electrolyte* fuel* and air temperature profiles were derived. 

For the three dimensional analysis of the stack* it was assumed 
,that all of . the walls except for the wall from tdiieh the air enters were 
maintained at a constant temperature. The rate of heat generation per 
.unit volume of the stack was assumed constant. An analytical solution 





I for the temperature profile was developed, assuming that the electrolyte 

s, 

and gas temperatures were not very different. 

Another paper (Ref. 15 ) considered various limiting and special 
cases to determine the maximum temperature of a stack. Two dimensional 
heat transfer analysis was carried out in the case of a thick stack 
where heat transfer in the direction of stacking was neglected. In the 
case of thin stacks, three-dimensional heat transfer was considered with 
each wall at a different temperature. Infinite series solutions were 
developed for both thick and thin stacks. The authors estimated the 

maximum stack temperature for the constant wall temperature case. An 
approximate formula to predict the effect of conductivities, size, and 
current density on the maximum stack temperature was developed. A 
( generalized analysis, which can incorporate the effect of finite 

y 

f- 

resistance to heat transfer at the wall, the effect of cold or hot 

i . 

feeds, or nonuniform heat generation, was also carried out using the 
method of Green's function. 

2. 2. 3. 2 Temperature Distribution 

The temperature distribution for the module was developed from the 
temperature distributions within representative slices or strips within 
a set of cell and cooling plate cells. The analysis includes conduction 
within bipolar plates, conduction between plates, the separate cooling 
effects of the process air and the coolant (basically air Is considered 
^ the coolant), and the temperature chax%e of air flows along their 
respective channels. The distribution of the heat generation- is 
determined from the current density distribution. 
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The model assumes that (1) the temperature gradients in the 
direction of the fuel flow are small. This assumption is Justified 
since the oiajor temperature gradients are in the air flow direction and 
since the heat capacity of the fuel stream is only a few percent of the 
heat capacity of the air stream. (2> The edge of the cell is operating 
adiabatically. (3) A half set of cell plates between cooling plates is 
analyzed, which includes one half cooling plate and two and a half cell 
plates. Thus because of the symmetry all of the stack behaves 
similarly. The geometry of a representative slice (Lx X Ly X Lz) 
through the stack is shown in Figure 2-17. 

Hathematical Formulation 

The material balances of the fuel and the oxident have been 
presented in section 2.2.2. There arc four exiergy balance equations for 
the cell plate, cooling plate, process air. and coolant. 
cell on process air side in air flow direction 


T . 3T 


- Kx 


x+t 


d T 


i 9 Tp 
- E ■ 

Pp ay 


Pooling plate in coolant direction 


- (V* -V) I (2-2-17) 


t* Ky 


T 


9y 


a T 

- + 2Kx — 

2 ax 


Cc a Tc 
3C+tV2 Pc 


(2-2-18) 


air side 
d Tn ^ Sp 

^ Cp ( ) 




(2-2-19) 





Symmetry for 
5 Cell Plate/ 
Cooling Plate 






coolan't side 


) Tc 

I y ~ me Cc 


( T-Tc ) 


( 2 - 2 - 20 ) 


Boundary conditions 

X = 0 ST/ox = 0 

y = 0 3T/0y = 0 

X = Lx 3T/0X = 0 

y = Ly 0T/6y = 0 

y = 0 Tp = Tp# inlet 

y = 0 Tc = Tc, inlet 

where 


symmetric condition 
adiabatic assumption 
symmetric condition 
adiabatic assumption 


m = mass flow rate, Kg/hr-channel 
c = heat capacity, J/Kg-K 

Ky = effective thermal conductivity of cell in flow direction. 




J/hr-m-K 

Kx = effective thermal conductivity of cell on stacking 
direction, J/hr-m-K 

t = thickness of cell including fuel and air channel, m 
xl = effective conduction distance from plate to upper cell 
plate, n 

x2 = effective conduction distance from plate to lower cell 
plate, m 

‘p = pitch of channel, m j - . a 

xl' = effective conduction distance from cooling plate to upper 
, ceil, plate, m _ _ r-; - a''. , ... ^ 




Lx,Ly = height and length of one slice, respectively, m 


OmGJI^AL.-P^ElS 
GF POC^ 





65 


v« = -AH/2F, V 

t' = thickness of cooling plate, m 
h s heat transfer coefficient, J/hr-n^-K 
S = perimeter of the channel, m 
Subscription 

p process air 
c = cooling air 

These simultaneous ordinary differential equations and 
corresponding boundary conditions were solved by the finite-difference 
method. The final difference equations are in next subsection. 

Finite-Difference Model 

The energy balance on an internal element j for bipolar 

plate i (2ii<Nl) can now be written as (see Figure 2-17 > 


( (Tpi.J = (V*-V) li,3 


( 2-2-2 T) 


The energy balance on an internal element j <2iJiN-l) of the 


cooling plate i=l can be written as 

wh*re 


+ (2 ■ ■ + 2 ) T1. a 

- 0+1 - (2 


A Y= 
Kx 


XI ^ 


( 2 - 2 - 22 ) 


i . 

1 f 
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The energy balance on interior element J of the symmetric 

plate i=Nl is 


Kv*t <5 Kx >. _ 

AY* X2 ^ N'l ,j 


. (JSCi— ) T + (2- 
^ AY* '' N1.0-1 

AY* ' m,j+1 - X2 ^ 




( 2 - 2 - 25 ) 


where S ~ 2 for odd values of NK 
S ^ 1 for even values of NK 

NK: the number of cell plates between two adjacent cooling plates 
Nl: 1 + NK/2 for even NX 
1 + (NK+D/2 for odd NX 


The energy balance on elements J=1 are obtained as above except 
for: the values of Ti,0 are replaced by Ti,l; the values of Tpi,0 are 

I 

replaced by TPO, which is the inlet process air temperature; and TCo is 
replaced by TCO, the inlet cooling air temperature. The energy balances 
on elements J=N are obtained from the above with Ti,J+l replaced by 
Ti,N. 


For the process air flow, one can set up NlxN equations of the form 
TPi,J = TPi, J-1 + (Ti,j - TTi,J) (l-e'^'’"''^) (2-2-2<l) 

Inhere 

'v h . \ Sp / ^ 

> » «-2-2S) 

® Ifti Cp - . 

For the cooling air flow, one obtains N equations of the form 



where 




Me Cc 


AY 


(2-2-27) 


Thus the total number of temperature equations matches the number 
of unknown temperatures and the set can be solved using the Gaussian 
elimination method with calculated or input values of cell voltages, 
current densities, mass flow, heat generation and heat transfer 
coefficients. Each resulting temperature distribution is used to 
recalculate the current density distribution until convergence is 
reached. The relationship between voltage and current and the 
calculation of heat generation have been presented in Section 2.2.1. 

Heat Transfer Coefficients 

An empirical equation (Ref. 18 ) for the Nusselt number for fully 

developed laminar flow in a rectai^ular channels is: 

Nuf = 3.61 + 4.63 (1-flO^-^ (2-2-28) 

where c4 = a/b; a is the smaller side of rectangular channel and 

b is the larger side of the channel. 

Near the inlet of a channel, the heat transfer coefficient is larger 
the fully developed value due to development of the laminar 
**®«ndary layer. If R is the ratio of the average Nusselt number for the 
region 0 to x to the fully developed Nusselt number, then (Ref . 19 ) 


?~-H » 1 I 

1 + 0.04 


(2-2-29) 
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where GZ: Graetz number = Re Pr 

Re: Reynolds number based on Dh 
P r: Prandtl number of gas 
Hydraulic diameter, n 

For turbulent flow, the average Nusselt number over the region 0 to x is 
described as (Ref . 22 ) 

Nu^ = 0.116 CRe^^ -125] Pr II + (D/x)^^^ ]. (2-2-30) 

The flow chart of the executive program for calculating the 
temperature distribution in the stack is shown in Figure 2-18. 

Results 

A temperature distribution in the fuel cell stack with 0.30 m x 
0.43 ffl cell plate has been determined using the computer program. This 
distribution is presented in Figure 2-19 as a field plot of lines of 
constant temperature on the second cell plate. The operating coxiditions 
are shown in the figure where Kx, Ky represent thermal conductivity 
along stacking direction and on the cell plate respectively. The 
maximum temperature occurs at the corner of fuel inlet and air outlet 
and the minimum occurs at the opposite corner. The stacking direction 
temperature profile will be discussed in the next section. 

Comparison with Experimental ‘Data ' . c rsr:rp 
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A comparison of the temperature profile on the cell plate as predicted 
hy the developed model was made with the Uestinghouse measured data 
(Ref. 11). Figure 2-20 shows the Uestinghouse measured temperature 
distribution with a ’’branch" cooling configuration which will be 
discussed in detail in the next section. Figure 2-21 shows the 
calculated results. The significant difference between the operating 
temperature and the anode inlet temperature in the experimental work 
conflict with the assumption which neglects the anode flow effect 
(Section 2. 2. 3. 2). Therefore, the computer peogram was modified to 
include a fuel flow cooling effect and the results are shown in Figure 
2-22 with the parameters listed on the figure. The ’’hot spot" shifted 
to right with the addition of the add anode flow effect (compare Figure 
2-21 with Figure 2-22 ). 


The temperature contour lines are very sensitive to many parameters 


which will be discussed in next section. For example. Figure 2-23 shows 
the contour lines for the same parameters as Figure 2-22 except the 


temperature of the process air (TKA) (not the cooling air) was decreased 
by 10®C. Therefore, considering the complex interaction among 
temperature levels, reaction rates, reactant depletions, current 
densities, heat generation rates and reactant distributions and that 
only temperatures can be measured directly; the agreement between 
calculated and measured temperatures is very good. 


2. 2. 3. 3 Parametric Sensitivity and Cooling Scheme 


The plate temperature is a function of the current density, the 



Westlnghouse experimental 







Calculated temperature dletribution - H2 not considered 







Calculated temperature distribution with 







The temperature distribution for tvo different sets of dimensions 
of cell plates were simulated and compared with, and the results were 
shown in Figure 2-19. All of the operating conditions were the same as 
in Figure 2-19 except for the plate dimensions (width and length). The 
first one (Figure 2-25 ) with the reversed width ant3 length of Figure 
2-19 and the second (Figure 2-26 ) with square dimensions. The results 
show that the rectangluar shape with a fuel flow length longer than the 
flow length will obtain the lowest peak temperature and have a good 
average temperature. Lower. peakTtemperatures ' mean a more uniform 
temperature distribution and a lower average temperature means that less 
auxiliary power will be required to pump the coolant. 





Temperature distribution 











Figure >2-26; Temperature distribution with square dimens ions (compared to Figure 2-19) 




Size 


The computer analysis was done for a sequence of cell plate sizes 


under the same operating conditions. The various sizes are 2” x 2" 


(5.08 cm X 5.08 cm), 6" x <15.2<1 cm x 10.16 cm), 10" x 7" (25.4 cm x 


17.8 cm), and 12" x 10' (30.48 cm x 25.4 cm). The normalized current 


densities and temperatures along the plate diagonals, are given in Table 


2-4 . The table shows that for a relatively small cell, both current 


density and temperature vary almost linearly from the reactants inlet 


corner to the exhaust corner. ks cell size Increases, both 


distributions become somewhat parabolic along the "k" diagonal, while on 


The "B" diagonal, the temperature extremes are accentuated for the 


larger plate sizes. 


The smallest plate shows the most uniform temperature distribution 


but the worst current density distribution. The greater change of 


hydrogen or oxygen concentration in the small plate results in a varying 


current density. 


The indicated trend in temperature extremes as plate size increases 


(while average temperature stays constant) can create a potential 


problem for very large plates. Some portions of the plates could 


possibly exceed the boiling point temperature of the electrolyte. 


It is noted that the _averagej^ current density and the average 
tea^rature are the same along the "k" and "B" diagonals (and the plate) 
for all of the plates sizes. This result indicates a convenience of the 




Diagonal Mode Ciirrent Densities on “A” Diagonal Current Densities on “B” Diagonal 

Coordinate Number (A/sq.cm.) (A/sq.cm, ) 

X/L 2^' X 2« 4« X 6« 7” X 10" 10" X 12" 2” x 2" 4” x 6” 7” x 10" 10” x 12" 
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Table 2-4 Effect of size on fuel cell performance 
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diagonal analysis method. 

Average Current Density 

Uhen the output power requirement of a PAFC system are determined, 
the average current density is also decided due to its relationship with 
voltage. The typical relationship between power and current density is 
shown in Figure 2-27 for a fixed plate area. In addition, examination 
of the effect of current density on the temperature distribution will 
decide the size or the number of plates needed. 

The computer program was run with four different values of the 
average current density which are 0.1, 0.2, 0.3> and 0.4 A/cn* . The 

diagonal analysis of "A" type diagonal (see Table 2-4 ) on a 25.4 cm x 
17.8 cm plate is - shown in Figure 2-28 and Figure 2-29 . It can be 

observed that the smaller average current density the more uniform the 
temperauture distribution. Another important result is that the current 
density distribution corresponded to the temperature profile for the 
larger average current density. This is contrary to the isothermal 
current density distribution which has a minimum value at the right side 
comer (Figure 2-13 ). 

Coolant Flow Rate 

A majority of the heat energy removal is accomplished by the flow 
of coolant through the cooling plate. For example;"-. 2 times the 
*toi^hiometrically required ^;.amount <2 stoieh) of c^hode gas is passed 

^hrwgh individual cells and a 30 stoieh catl^e gas is going through 

52 'j' 
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e 2-29 Effect of the mean C.D. on the diagonal temperature distribution 
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cooling plate. From a cost view point, the greater the flow rate of 
coolant the greater the auxiliary power needed to recycle the coolant. 
Horeover, the flow rate of coolant will also affect the temperature 
profile. Figure 2-30 shows an example for the stacking direction 
temperature profile with two different stoich values of coolant. The 
more coolant that is used, the lower the mean temperature will be, but 
the greater the temperature difference between plates will be. 


The thermal conductivity in flow direction and inlet temperature of 
both fuel and process air also affect the temperature distribution. 
This will be discussed after we compare different cooling schemes. 


Cooling Scheme 


The cooling factor is conceived to be a function of the heat 
transfer characteristics, plate size, and stack construction. The 
latter will primarily specify the number of power plates between a pair 
of cooling plates. The heat transfer characteristics will be a funxrtion 
of the type of coolant (gas or liquid), cooling plate design, and the 
thermal conductivities of the plate material. Some of these factors 
have been examined in the previous subsection. This subsection will 
concentrate on the design of the cooling chaxmel. 


There are three configurations of cooling channels considered here, 
**ose nomenclature and definitions are ais follows, 

1. Straight: the dimei^ibns of cooling channel .are :fixedi 

2. Branch: the cooling channel is branched along the coolant flow 

direction, one example is in Figure 2-31 . 


Plate Number Along Stacking Direction 

2-30 Con 5 >arison of the effect of the coolant flowrate on the 
mean temperature for plates in the stack direction 
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3. Varying width: the width of cooling channel is different along 

the fuel flow direction. 


Hfter the cooling stream has become fully developed the heat 
transfer coefficient drop dramatically (see Subsection 2. 2. 3. 2). The 
"branch'’ configuration was designed to prevent the formation of fully 
developed flow and to increase the flow rate (as the total crossectional 
area is decreased). The "varying width" configuration will put more 
coolant on the larger heat generation side, but the heat transfer 
coefficient does not change. 




In the "branch” configuration, a study was done to optimize the 
arrangement of different sections. The number of sections and the width 
of each section were hold fixed. The computer code developed in the 
previous section was modified to meet this configuration. Table 2-5 
presents the results for various length ratios (dimensionless) with 
comparison to the "straight" configuration. The estimated standard 
deviation of the temperature and the peak temperature provide an 
Indication of the uniformity of the thermal distribution. The average 
temperature can be used to estimate the cooling effect. Both of these 
results are not quite concordant. For example. Figure 2-32 plots the 
results for different length ratios of the first section, which include 
the "straight" configuration. It csm be observed that the most uniform 
distribution occurs at higher (around 7) ratios, trtiereas, the best 
cooling effect occurs at lower (around 2) ratios. The results of the 
comparison among different section length arrangements are shown in 
.^ifWre 2-33 in which the temperature profiles along the coolant 


n 


lli^ : 



First Separating Section Eatio ( 10 means no separation;) 

Figure 2-32 Con 5 )arison of ten 5 )erature on the cell plate with 
different first separating section • 
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direction are drawn. There are five arrangements in this figure, they 
are the straight configuration, the longest first section (7:2:1), the 
longest second section (1:7:2), the average length (4:3:3), and the 
longest third section (1:2:7). The arrangement of average length shows 
the best cooling effect and good thermal uniformity (the values are 
listed in Table 2-5 ) . hore detailed comparisons among the three 
configuration are shown in Table 2-6 with a fixed ratio of coolant and 
Table 2-7 with a fixed operating temperature. The section ratio in 
"varying width" configuration is the ratio of the number of channels in 
each section which have different width ratios. The different channel 
width will result in different coolant flow rates. 


It can be observed that: (1) from the uniformity view point, the 
branched channel is 23% ( (10.-7. 7)/10. ) better than "straight" 
configuration and 145J ( (8.9-7.7)/7.7) better than "varying width" 
configuration and the peak temperature is much lower (Table 2-6 ); (2) 
from the cooling effect view point, the branched channel is 17% 
((35-29)/29) better than the other two configurations (Table 2-7 ). 


u 
0 ! 
i ? 


The "varying width" configuration requires further examination of 
the effects of the section and width ratios on the uniformity to 
determine an optimal design. However, there is no optimal ratio for all 
operating conditions. For example, the results for the calculated 
uniformity, compared to the branched channel are quite different in Table 
2-6 and Table 2-7 . Although there Is no improvement in the cooling 
effect, the "varying width" type is attractive due to its simplicity for 
manufacture and reliability over the "branch" configuration. 



no separation 

branched channel 
(section ratio 2:3:5) 

varying width 
(section ratio 3:2:1) 
(width ratio 1:0.95:0.9) 

^•isated 

i*j«^ard 

jeriation 

0,0236 

0.0202 

0.0152 

»T«rage 

usperature 

180.0 

171.1 

180.2 

}*Ak 

u^erature 

17.1 

8.0 

13.3 

t#s? era tore 
tatiaated 
ttandard 
ieriation 

10.0 

7.7 

8.9 


tempera-bjre; ®C 
C.D.; An^j/cm^ 


C«n«ral input data: 

. 2 
J.Ate dimension: 12, x 10, in 

utilization of H 2 : 0,75 

utilization of O 2 : 0,50 

Ifiput H2 fraction: 0,75 

icput cooling air ten^jerature: 575 

4y«ra^ C,D,: 0,300 Amp/cm^ 


Table 2-6 Comparison of different cooling configurations 
with fixed ratio of cooling air 
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no separation 

branched channel 
(section ratio 5:2:1) 

varying width 
(section ratio 3:2:1) 
(width ratio 1:0,91:0.77) 

:,o. 

♦itimated 

standard 

jeviation 

O.O8I72 

O.O615I 

0.05382 

iverage 

terperature 

190.0 

190.2 

190.8 

;«aic 

•.•rperature 

36.3 

24.7 

25.6 

temperature 

eatiaated 

standard 

‘rviation 

18.77 

13.52 

15.62 

ratio of 
tooling air 

35 

29. 

35 


sniti teniperature: ®C 
C.D.: An^j/cm^ 


funeral input data: 

plate dimension: 15.95 ^ 12. in^ 

-tilization of H2: 0.8 

-tilization of O2: O.5 

i2jmt H2 ratio; O.76 

Jap^t CO ratio: 0,01 

iJ5put cooling air ten^serature: 405* 3 *2 

•^erage C.D.: O.325 Amp/cm^ 


Table 2-7 Comparison of different cooling configurations 
with the Same average operating temperature 



Thermal Conductivity in The Flow Direction 


the temperature profiles along the coolant flow direction are 
simulated in Figure 2-34 for two different Ky values (thermal 
conductivity in flow direction) in the "branch” configuration. For the 
smaller Ky value, the break of the fully developed phenomenon can be 
observed more apparently which increases the heat transfer coefficient. 
Ks was expected, the greater Ky values produces better thermal 
distribution and a worse cooling effect. The material and construction 
of cell plate decide the Ky value. 

Inlet Temperature of Process Mr 

Though the cooling effect of the process air is less than 1C% of 
the total cooling, the inlet temperature of process air (TKi^), or more 
precisely the temperature difference (AT) between average operating 
temperature and TKk, is an important factor in thermal analysis. For 
example, in the "branch" configuration, a small AT (183- 155=28 "O will 
obtain more uniform temperature profile than that of higher aT 
(163-118=45 *0. The results are shown in Figure 2-35 and Figure 2-36 . 
Therefore, it is suggested that both the anode and cathode inlet 
temperatures should be close to the average stack temperature. 


I 


s 


I 



Figure -2-54 Comparison of distribution with different 
thermal conductivities in the plate direction 


Ky* 50 
Ky= 15 


9 ' 



Teaperature 
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2.3 Steady State Performance of PAFC System 


The computer codes developed in previous sections were combined to 
simulate the PAFC system performance. The main program is the executive 
which controls the flows, manages and tests the iteration procedures, 
reads input data, and writes the final results. There are two PAFC 
systems which will be simulated: The CSU design and the Uestinghouse 
design. 


Results for the CSU Design 


The lumped model of each component was used to simulate the 

original CSU design (Figure 2-1) except in the reformer where the 

distributed model was used. The flowchart of the executive program is 

shown in Figure 2-37 and the input data definitions and values are 

listed in Appendix 5. Figure 2-38 provides the final results which are 

presented as the P-T-V <V as a flow rate) status of each stream numbered 

on the flow diagram. The calculated electrical and heat generated are 

0 

12^ Kw DC and 4.78 x 10 J/hr, respectively with a fuel cell efficiency 
of 37.45J. 


There has been much interest in the effect of alternate commercial 
fuels on the performance and costs of PAFC power plants. The computer 
program developed can be modified to allow for methanol or naphtha as 
li^iut fuel. The system with methanol input fuel obtains the highest 
effeciency (approximately 42% for total the system) among the three 
fuels. IJore detailckl discussions on this topic are presented in Ref. 


divider 
CALL LIVID 




2—37 Plow chaurt of executive program for simulating CSU's PAPC system 
steady state performance 
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Wgure 2-^38 Temp era ttire, preaeure, and mass balance of Cb'U’s PAJ*X! system 


105 


Comoerison with the Uestinghouse Conceptual Design 

Both the lumped and distributed models were used to simulate the 
Uestinghouse PAFC system (Figure 1-3). The estimated dimensions of the 
heat exchangers, shift converters, and reformer discussed in previous 
sections were used as input data for the distributed model simulation. 
The definitions and values of input data are listed in Appendix 6. 


>■* * 
m- 

Si 



Because the flow diagram of this system differed from the CSU 
system, the executive program was modified to meet this change. The 
flowchart of executive program is shown in Figure 2-39 . Results for 
the mass balance, temperature, and pressure of the gas stream at various 
locations for Uestinghouse PAFC system are given in Table 2-8 . The 
stream number refers to Figure 1-4 and the gas component values are in 
mole fractions. There are three values presented for each stream, the 
first row is the Uestinghouse conceptual design (Ref. 12), the second 
row is the results from the "lumped" simulation, and the third row is 
the results from the "distributed" simulation. A comparison of results 
is shown in Table 2-9 . Because the type of reformer in the distributed 
simulation is "once-through" and the reformer in the Uestinghouse 
conceptual design is "regenerative"*, the outlet temperature of 
reforming gas is higher (939*K) and the temperature of combustion gas is 
lower (835"K) in the CSU simulation than that~ (867*K and lOBl'K 
respectively) in Uestinghouse conceptual design. The results in both 
shift converters are also affected by this temperature difference. One 
serious problem that arises with the (SU "once-through" reformer is that 
not enough heat for evaporating water is available from the outlet 





/READ 

input data/ 


calculate 
conditions of 
direct cOTitacI 
condenser 


asstune 
X and DTM2 


I calculate 
concentration 
of anode 
inlet 


calcxilate 
input flow ra ;e 
of Stack 
CALL DATACB 


calculate 
input fuel 
rate 

fuel input 
con5)ressor 
CALL COMP 



^gure 2-39 The flow chart of executive program for steady state simulation 
of Westinghouse PAEC system 
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CO 

CC2 

U20 

H2 

H2 02 

rUM RATt(lb/hr) 

Ttre>HIATUHi;( F) 
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0 

0 
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.010 
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14084.68 

80.00 

14.70 

0 

. 0 
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.010 
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.782 .208 

14326.57 

80.00 

14.70 

0 

0 
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.010 
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.010 
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14084.68 
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14326.57 
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0 
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.010 
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0 
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.010 
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14084.68 

214.45 
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.010 
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14326.57 
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18.80 
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.782 .208 

14358.00 
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0 
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.782 .208 

14084.68 

747.61 
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0 
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.010 
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.782 .208 

14526.57 
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0 

0 

0 

.010 

0 

.782 .208 

14358.00 

700.00 

20.00 

.950 

0 

0 

0 

0 

.050 0 

2855.75 

80.00 

49.98 

.955 

0 

0 

0 

0 

.047 0 

2800.00 

80.00 

49.98 

.950 

0 

0 

0 

0 

.050 0 

2797.00 

77.00 

50.00 

.950 

0 

0 

0 

0 

.050 0 

2853.75 

181.42 

98.48 

.953 

0 

0 

0 

0 

.047 0 
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183.83 

99.96 
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0 

0 

0 

0 

.050 0 

2797.00 
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0 

0 

0 
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0 
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0 
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700.00 
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.950 

0 

0 

0 

0 

.050 0 

82.64 
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0 

0 

0 

0 

.047 0 

81.09 

724.72 
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.950 

0 

0 

0 

0 

.050 0 

81.00 
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.281 

0 
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.704 

0 

.015 0 

9907.54 

553.36 

97.80 

.282 

0 
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.704 

0 

.014 0 

9754.09 

551.33 

1C2.;-- 

.271 

0 

0 

.714 

0 

.014 0 

9779.00 

556.00 

100.00 

.281 

0 

0 

.704 

0 

.015 0 

9907.54 

856.50 

97.88 

.282 

0 

0 

.704 

0 

.014 0 

9754.09 

798.80 

102.26 

.271 

0 

0 

.714 

0 

.014 0 

9779.00 

800.00 

100.00 

.014 

.101 

.070 

.222 

.583 

.010 0 

9907.55 

1229.84 

76.04 

.012 

.122 

.051 

.237 

.569 

.009 0 

9754.09 

1100.06 

74.97 

.014 

.119 

.053 

.242 

.562 

.009 0 

9779.00 

1100.00 

75.00 

.014 

.101 

.070 

.222 

.583 

.010 0 

9907.55 

990.57 

76.04 

.012 

.122 

.051 

.237 

.569 

.009 0 

9754.09 

899.19 

74.97 

.014 

.119 

.055 

.242 

.562 

,009 0 

9779.00 

899.00 

75.00 

.014 

.101 

.070 

.222 

.583 

.010 0 

9907.54 

792.70 

75.84 

.012 

.122 

.051 

.237 

.569 

.009 0 

9754.09 

721.99 

74.96 

.014 

.119 

.053 

.242 

.562 

.009 0 

9779.00 

720.00 

75.00 

.014 

.063 

.108 

.164 

.621 

.010 0 

9907.53 

887.46 

71.71 

.012 

.060 

.113 

.175 

.632 

.009 0 

9754.09 

651 U4 

74.96 

.014 

.058 

.115 

.180 

.623 

.009 0 

9779.00 

840.00 

70.00 

.014 

.063 

.108 

.184 

.621 

.010 0 

9907.53 

411.87 

71.71 

.012 

.060 

.113 

.175 

.632 

.009 0 

9754.09 

442.41 

74.96 

.014 

.058 

.115 

.180 

.623 

.009 0 

9779.00 

397.00 

65.00 


Table 2-8 Ten 5 )erature, pressure, and mass balance of Westinghouse PAPC system 
with three kinds of simulation 


n 

.014 

.002 

.169 

.123 

.603 

.010 

0 

9907.52 

557.25 

66.05 


.012 

.009 

.164 

.124 

.605 

.009 

0 

9754.09 

553.65 

74.95 


.014 

.013 

.160 

.155 

.669 

.009 

0 

9779.00 

495.00 

65.00 ■ 

18 

.014 

.002 

.169 

.125 

.683 

.010 

0 

9907.51 

481.23 

66.02 


- .012 

.009 

.164 

.124 

.683 

.009 

0 

9754.09 

491.73 

74.95 


.014 

.015 

.160 

.135 

.668 

.009 

0 

9779.00 

401.00 

60.0c 

19 

.014 

.002 

.169 

.123 

.663 

.010 

0 

9907.51 

191.31 

66.02 


.012 

.009 

.164 

.124 

.603 

.009 

0 

9754.09 

211.02 

74.95 


.014 

.013 

.160 

.135 

.669 

.009 

0 

9779.00 

200.00 

60.00 

20 

.016 

.002 

.168 

.024 

.759 

.011 

0 

8353.79 

123.20 

60.01 


.012 

.009 

.167 

.063 

.741 

.009 

0 

8879.36 

159.20 

60.01 


.015 

.014 

.171 

.077 

.715 

.010 

0 

8826.00 

159.00 

60.00 

21 

.016 

.002 

.188 

.024 

.759 

.011 

0 

8353.79 

450.72 

60.01 


.012 

.009 

.167 

.063 

.741 

.009 

0 

8879.56 

449.27 

60.00 


.015 

.014 

.171 

.077 

.713 

.010 

0 

8826.00 

576.00 

60.00 

22 

0 

0 

0 

1.00 

0 

0 

0 

7136.45 

525.21 

95.72 


0 

0 

0 

1.00 

0 

0 

0 

7035.18 

328.56 

100.45 

23 

0 

0 

0 

1.00 

0 

0 

0 

7136.45 

441.98 

95.52 


0 

0 

0 

1.00 

0 

0 

0 

7035.18 

445.77 

99.95 


0 

0 

0 

1.00 

0 

0 

0 

7063.00 

480.00 

100.00 

24 

.040 

.006 

.479 

.062 

.386 

.028 

0 

7418.01 

845.90 

47.48 


.029 

.022 

.410 

.153 

.564 

.023 

0 

7891.28 

760.12 

49.97 


.035 

.051 

.398 

.180 

.532 

.025 

0 

7924.00 

. 761.00 

50.00 

25 

0 

0 

.221 

.235 

0 

.531 

.015 

21585.54 

5266.50 

25.49 


0 

0 

.207 

.269 

0 

.512 

.012 

22296.94 

3176.50 

22.91 


0 

0 

.207 

.270 

0 

.511 

.012 

22564.00 

5184.00 

15.00 

26 

0 

0 

.221 

.235 

0 

.551 

.015 

21585.34 

1042.89 

21.14 


0 

0 

.207 

.269 

0 

.512 

.012 

22296.94 

1279.95 

20.62 


0 

0 

.207 

.270 

0 

.511 

.012 

22364.00 

1432.00 

15.00 

27 

0 

0 

.221 

.235 

0 

.551 

.015 

21585.34 

525.21 

20.51 


0 

0 

.207 

.269 

0 

.512 

.012 

22298.94 

519.56 

20.00 


0 

0 

.207 

.270 

0 

.511 

.012 

22364.00 

446.00 

15.00 

26 

0 

0 

.221- 

.235 

0 

.531 

.013 

21585.34 

249.92 

20.41 


0 

0 

.207 

.269 

0 

.512 

.012 

22298.94 

231.09 

19.90 


0 

0 

.207 

.270 

0 

.511 

.012 

22364.00 

555.00 

15.00 

29 

.950 

0 

0 

0 

0 

.050 

0 

2771.11 

735.48 

95.53 


.953 

0 

0 

0 

0 

.047 

0 

2718. 91 

724.72 

99.95 
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Table 2-8 continued 
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Westinghouse 

Our 

Oir 


Conceptual 

Oeai0i 

Besults-1 

Besulte-2 

Systao SiBulatlon Method 


'Rieraodynaaic 

: Distributed 

Katural Cos Comp. Vol % 

-CH4 

90 

95.3 

95 

-C2H6 

5 

0 

0 

-N2 

5 

4.7 

5 

•Preeeurc peim 

100 

100 

98.5 

Bafoner Tuba ^rpa Baganarmtlva 

Bagenerative 

Once- through 

-Approach Slit, Temp. *p 

25 

25 


•Caivarslon of NathaDe 

0.925 

0.935 

0.924 

-Baformiisg Gas Inlet Tea^. *P 

300 

798 

856 

-aefoming Gaa Oitlat Temp. *F 

1100 

1100 

1230 

-Heforming Caa Inlet Prasaura psia 

100 

102 

98 

-Beformlng Gaa Oitlat Preaaura pala 

75 

75 

76 

-Combuation Gaa Inlet Tamp. *P 

3184 

3176 

5267 

-CombuBtion Gas Outlet Temp. *F 

1432 

1280 

1043 

High Temp. Shift Converter 

-Approach Slit, 7mp, *F 

-25 

-25 


-Converaion of CO 

0.51 

0.5 1 

0.38 

-Inlet Tamp. *r 

T20 

"22 

793 

-Outlet Temp. *F 

840 

851 

336 

-Inlet Preeaure pala 

75 

75 

76 

-Outlet Preasure pala 

70 

75 

72 

Lay Temp. Shift ConvarteT 

-Approach Oiff. Tamp. *F 

-25 

-25 


-Converaion of CO 

0.776 

0.35 

C.?7 

-Inlet Temp. 'F 

397 

442 

412 

-<Xitlet Temp, 'r' 

495 

554 

557 

-Inlet Preasure paia 

65 

75 

72 

-Outlet Preasure psia 
Fuel Cell Stack 

65 

75 

<6 

-9eai^ Currant Density mA/cm^ 

325 

325 

325 

-Deal^ (Calculate) Cell Voltage 

0.683 

0.616 

0.665 

volts/cell 

-Operating Tei^erature (Average) 'T 

376 

376 

375 

-Operating Pressure pala 

50 

50 

50 

-Effective Call Area em^ 

1080 

1080 

103c 

-Anode Inlet Tengierature ‘? 

376 

449 

451 

-Anode Outlet Tea^arature *? 

376 

575 

562 

-Cathode Inlet Tamperatore *F 

3- 

368 

350 

-Cathode Outlet Temperature *F 

376 

375 

395 


Table 2-9 Comparison of simulations of Westinghouse PAPC system 
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cosmustion gas at This condition lacks 1/3 to 1/41 the heat duty 
(see the output for this simulation in {Appendix 7). One way to amend 
this problem is to recover some heat from the product gas by adding a 
heat exchanger between the evaporating water and product gas. 

(isnsidering the complex interactions of temperature and composition 
at each location, the agreement between the CSU simulations and 
Westinghouse conceptual design is very good. 


•I’egenerative ref brmer';- • the 'j^bduct ^ gas is reveled through ^^the ■mbst 
inner tube to" recover somd energy (heat).’ This design • will’ save icr-i5% 
®^ "the heat lieededi in ref ormihg (Ref. 13) ' -■ 

h:? 1 r .r- ^ ~ vs* ... .... 
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2.^ Process Synthesis with Keat Exchanger Network 


The CSU PRFC system developed originally can be cost improved by 
considering the energy conservation. Any heat recovered and reused in 
the process not only reduces the amount of fuel consumed but also 
reduces the amount of heat rejected to the cooling system. One way to 
solve this kind of problem Is to find the optimal heat exchanger 
network. This task consists of finding a feasible sequence of heat 
exchangers in which pairs of streams are matched, such that the network 
is optimized as the total cost is minimized. 


A new combined algorithmic-heuristic approach has been developed. 
This three-phase algorithm can generate, with great ease and 
considerable speed, a near optimal exchanger network. The description 
of this new algorithm and a comparison with other approaches is 
presented in Appendix 3. This approach was used to improve the heat 
exchanger network in the CSU original PAFC system shown in Figure 2-1. 
Figure 2-<IO shows the network with the corresponding temperatures. 
Here the cooling media (steam without phase change) in the fuel cell 
stack is to be cooled from ^33®K to <J13®K by evaporating the pressurized 
steam (steam no. 25) and the outside loop cooling water. The outlet 
cooling water can be used in domestic applicitions. Figure 2-41 shows 
the improved heat exchanger network for the same operating conditions. 
The cost comparison is shown in Table 2-10, where the assumptions for 
the heat exchanger cost calculations are the same as in the sample 
problems described in Appendix 3 Table 2. There is 10% cost improvement 
Tor the new design over the old one. A greater reduction in cost could 
expected if the inlet temperatures of the reformer, shift converters. 


SHIFT 



Figure 2-40 Heat exchanger ne^ork of original PAFC system 






ten5>erature : K 



Figure 2-41 Optimal heat exchanger network of PAPC system 






Cost of 


Original network 
(Figure 2-40 ) 


Optimal network 
(Figure 2-41 ) 


E - 1 
E - 2 
E - 5 
E - 4 
E - 5 

E - 6 

E - 7 

E - 8 

E - 9 
E - 10 
C - 1 
C - 2 


Total 

Table 


83. 

251. 

164.8 

149.5 

175.1 

66.4 

529.6 

209.1 

1216.5 (cooler) 

112.6 

176.7 

92. 

132. 

191.2 

972.8 

507.1 

115.6 

2071.7 (cooler) 

2779.5 

411.5 


5522.1 


4885.6 


;-10 Cost comparison of heat exchanger network 



and mixer were no* fixed. 


The resul^s presented depend upon the assumptions made 
estimation and heat capacity calculations. However, this 
provides an optimal scheme for energy conservation. 


in cost 
approach 



CHAPTER 3 


SinULRTICi'J OF TRAJ'JSIENT STATE WITH LOAD CHAJ.’GE 

In normal operation, the plant control system will act to maintain 
the electric power output at a level required to meet the plant load 
demands. The cathode air and anode fuel gas flow rates will be 
maintained at prescribed levels corresponding to the current demand 
level. The cooling air flow will be adjusted to maintain a constant 
fuel cell temperature. The fuel cell pressure level will be adjusted to 
match the compressor characteristics. 

The daily electrical energy consumption either at a commercial or 
residential site is varying. Load change is an important and frequent 
operation in the power plant. Since the PAFC system can be subjected to 
sudden load changes and load ramping, an understanding of the effects of 
these transient conditions on the PAFC system’s performance is essential 
for the optimal design and control of the system. For the PAFC system, 
the controlled variable is the level of output AC power required, and 
the manipulated variables are the input fuel and air quantities. In 
addition to the orderly analyses, evaluation of the performance under 
catastorphic conditions can provide a means of designing emergency 
shutdown systems, and thus prevent disasters in PAFC plant. Basically, 
the problem is a "what would happen if" situation that must be analyzed, 
P*^®Perably, before the plant is built. 
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Therefore, the pro=e« dynaeics Ihvolve the hodelins of the 
.rooessln, tytteh Ih order to optlrdte thvestoeht and operatihr costs. 
,,Wev. autoaatio control, and provide for startup, shutdovn and 

emergency procedures. 
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3.1 Modeling Procedures of the Transient Response of the Syster, 

The steps to model the transient response of the PAFC syster. are 
described as follows: 

(1) Set up the mathematical model for each major fuel cell system 
component . 

(2) Estimate the design dimensions for each major fuel cell system 
component . 

(3) Simulate the distributed system performance at steady state as the 
initial conditions for the next step. 

«1) Calculate the transient state system responses. 



The results of steps (2) and (3) have been shown and discussed in 
Chapter 2. The following chapters will concentrate on steps (1) and 
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3.2 General Procedures for Transient State Simulation of Ka jcr 

Components 

For computer simulation of the dynamic system, a mathematical model 
must first be developed which defines the input-output relationships of 
a process as a function of time. Usually these mathematical 
relationships are in the form of simultaneous differential equations. 
These differential equations are determined either empirically, from the 
operating system or a similar operating system, or analytically, from 
the physical laws governing the dynamics of the system. The computer 
simulation described in this work was derived from a combination of 
these two techniques. 

A system of simultaneous partial differential equations was derived 
and subsequently solved from the physical and chemical laws, such as the 
law of conservation of mass, energy, aiid momentum, and the defining 
equations governing the PAFC system. 

In solving these types of problems, the first step is problem 
definition. Once the problem has been defined, the next step for 
solving it is usually the application of conservation equations. These 
apply to mass, energy and momentum balances. Generally, this procedure 
is applied to steady- state or static problems. The general equations 
for mass, energy or momentum conservation for a system are given by: 

■ 'Flow Flow i_'Re'actibh ’_ ■ v .. '..c.- 

In out ' product “ Accumulation , , (3.1) 

Energy^ Energy , Reaction _ AccuMlation 
' in ■” out ' ■ - eher^ ■ ~ ’ of energy 


(3.2) 
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Sum of forces acting = Change of momentum (3.3) 

In Eq. (3.1), the flow and accumulation terms can represent the 
amount of a particular component, or all (i.e., total) components of the 
system. In Eq. (3.3), the change of momentum represents the mass 
multiplied by its velocity. Eq. (3-3) is applicable to a specific 
direction-i.e. , it is a vector equation. 

The general procedures for transient state simulation of major 
components are as follows: 

(1) Set up the mathematical model for the transient condition. 

(2) Formulate explicit differential-difference working equations* from 
the mathematical model. 

(3) Delete time varying terms to obtain steady state solutions as the 
initial conditions. 

(q) At t=0, perturb the controlled variables with a desired change. 

(5) Solve these equations with a numerical integration routine for each 
time step. 

(6) Increase time by one time step, if steady state is reached plot the 
results; otherwise go to (5). 


•Differential-difference equation 

By substituting the distance varying terms with backward or forward 
differences and keeping the equation in explicit form, the equations 
thus obtained form a set of initial value ordinary differential 
equations which can be solved by standard integration methods. 





3.3 Concuter Solutions for Differential Equations 


The most powerful tool for solving the differential equations for 
modeling process-dynainics problems is probably the digital computer. 
Digital computers eliminate many of the disadvantages of their 
analog/hybrid counterparts, for example: 

(1) Initial setup time is less. 

(2) Once a digital program is developed that adequately simulates a 
certain phenomenon, we can retrieve the program and simulate the 
phenomenon under different conditions at a later time with less 
difficulty. 

(3) Eeproducibility of the results is nujch greater with the digital 


computer. 

(4) The digital computer can handle very small numbers and very large 
numbers. 


Problems arise in solving systems involing differential equations 
when numerical integration for digital simulation is used. All of the 
numerical integration methods that will be described require that a 
problem have known values of the dependent variables at some time, 0 , 


that is usually set equal to zero. He then predict the values of the 
dependent variables that satisfy our differential equation at successive 
time intervals via one of several numerical techniques. 


- t,-. ^ :.r r - retie:-: p-:r:.v 

Host generalized numerical-integration subroutines are set up to 
*olve equations of the following type (written for n first-order 


«*!uations): 




d9 = fl (0, yl, y2, yn) 
dy2 

d0 = f2 (0, yl, y2,..., yn) 


dyn . 

d9 = fn <0, yl, y2,..., yn) 


( 3 .^) 


(1) Euler method. He can calculate the point y at time (0+h) by Taylor 
series for Equation (3. 4) as: 


dy <0 ) h* d*y (0 ) 

y (0,+h) = y <0. ) + h d0 +2 d0* + ... 

where h is the step size. 


(3.5) 


In using this equation, the value of y(0. ) is given by the initial 
condition and y' (0, ) is evaluated from f (0. , y, ), given by the 
differential equation, dy/d0 = f <0. , y). It will of course be necessary 
to use this niethod iteratively, advancing the solution to 0 =0o + 2h 
after y (0,+ h) has been found, then to 0 = 0, + 3h, etc. 


Adopting a subscript notation for the successive y-values and 
representing the error by the order relation, we may write the algorithm 


for the Euler method as: 

yn+1 = yn +hy'n + 0(h’ ) 


(3.6) 


(2) Fomrth-order Runge-Kutta method. The most common used set of values 
in fourth-order Runge-Kutta method is summarized in Table 3.1. 

Runge-Kutta-Merson method. For inte^ation purposes, procedures 
vith automatic step-size control are recommended. The 
^®ge-Kutta-Merson routine is summarized in Table 3.2. The error of 
^***«gration can be (after oi»e step) estimated from the relation for e. 
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Provided E is the maxirjan permissible error, the step size is halved, if 
E<e, and doubled, if e<E/32. 


Table 3..I ?ourth-Crder ftunge-Xutta Soutine 

ynfl » + g ( lci + 25 c2 + 2 IC 3 ♦ ) + O(h^) 

» hf ( 6n, yn ) 

k2 • jh, yjj jki ) 

kj * hf ( On * jh, ♦ jkj ) 

:<4 » hf ( Gj. ♦ h, Vj. -► kj ) 




Table 3.2 The Hunse-iCutta-Kerson Soutine 

9 

^nfl - yn 2 ( * ^4 * ''S ) * 

■‘I - 3^ ( 9a. yn ) 

ic 2 • 3 ^ < 9j, ♦ y„ ♦ Ic, ) 

icj - ^f ( V y„ ^ |2 ) 

■ ■ 1=4 * 3^ ^ y„ ^ ^ ^ ) - 

^ -i =5 - 31 ( 9n + h, y„^* - fk* + 6 IC 4 ) 

J; Se' » ki - ^5 4^4 - |5 ; j- 


■this chsptieT 1.^ 


•KnCv:T: 5 n i' i ' * ■ ; a •*'. - 1‘ . 


CHAPTER H 


DYriAMIC SIHULATION OF THE FUEL PROCESSING SUBSYSTEIl 

The fuel processing subsystem, discussed in Section 2.1 for design 
under steady state conditions, vas analyzed dynamically for simulation 
of the transient response due to changing loads. This chapter describes 
the component and subsystem dynamic analyses. 

The subsystem considered here is basically from the Mestingbouse 
conceptual design (Figure "the analysis. More attention was 
placed on the process components. Therefore, the auxiliary components 
were assumed to be at steady state conditions and were neglected in the 
simulation of the transient state response. Figure 4-1 shows this 
simplified subsystem. 

h convenient method for analyzing the dynamic characteristics of 
dynamic in a component or system is to develop a computer program for 
the mathematical model and solve it. The objective here is to set up a 
system of simultaneous differential -difference equations, and use the 
numerical Integration methods discussed in previous chapter to solve 
-these equations. The general flowchart of the executive program used in 
this chapter is shown In Figure 4-2. ^ 

















record fjow 
rate of 
tube side 


record flow 
rate of 
shell side 


calculate time varying 
terms 

CALL FUl^CTH 

(heat exchanger) 
FONCTS 

(shift converter) 

FDUCTR 

(reformer) 


“^ 1 ( 1 ) ^ 


WRITE 

time 

varying 

term 


integration 

CALL 

EULER 


Wgur. 4-2 General flowchart for calculating transient change in heat exchanger, 
shift converter, or reformer 


P^E tS 

0 
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/g.l Keat Exchangers 

The heat exchangers in the PAFC system are simulated as counter 
current double-pipe heat exchangers as described in Chapter 2. The 
unsteady-state operation is considered with the dependent temperature 
variables as functions of the independent variables, time and distance. 
Since there is more than one independent variable, the relationship 
between temperature, time, and distance can be stated as a partial 
differential equation. 

^.1.1 Rathe matical Rodel 

The heat balance equations on the tube-side and shell-side can be 
written as follows: 

rate of heat accumulation * heat flow in - beat flow out + heat 

transfered 

tube-side: 

XPi(di/ Cp dt/ae = -fi vi^di* Cp St/^z + 4 hiXdl (Tw-t) 

M-I-1) 

shell-side: 

P© Ao Cpo 3l/d8 = -Po Ao Cpo vo dJ/dZ - ho do <T-Tw) «J-l-2> 

*diere the variables are the same as those used in the distributed model 
(for steady state conditions). t 

The equation for the temperature of the tube-wall separating the 
*hell and tube sections is 
tube-wall: 

pw<do*-dl*)Cpw dTw/d0 s ho do <T-Tw) - hi di <Tw-t) (4-1-3) 
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where the initial conditions 

Boundary conditions 

dt inlet _ Q 
d0 


are the steady state solutions 


dT inlet 


d8 


= 0 


(^- 1 - 4 ) 


In Equations (<l-l-l) and (4J-1-2) the inside and the outside 
composite coefficients include the thermal resistance of the tube-wall. 
If we add half of this resistance to both the inside and the outside 
heat transfer coefficients^ we can calculate the composite coefficients 
which are given by: (Ref. 21) 

hi = [l/hi' + R + di (do-di) ✓ kw (do+di)] 


ho s [l/ho' R + do (do-di) / kw (do+di )j * 
where hi' : inside tube heat transfer coefficient, J/m*-s-K 
ho' : outside tube heat transfer coefficient, J/n*-s-K 
R : fouling factor, m*-s-K/J 

kw : thermal conductivity of tubewall, J/m-s-K 
The fouling factor R was given the value 0.000176 m’-s-K/J. 

Differential Difference Eauations 


(4-1-5) 

(4-1-6) 


Since the objective was to reduce the mathematical model to a set 
of simultaneous ordinary differential equations (i.e. only one 
independent variable, time), the equations were *^finite difference" to 
a set of differential-difference equations. The heat exchanger was 
divided into tl sections: with each section on the tube-side corresponding 
“to a section on the shell side; a typical section J appears in Figure 
. Fluid on the shell-side flows from section J+1 to section J, 
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on the tube-side, it flows from section J to section j+1. The 
final set of differential-difference equations for the heat exchanger 
are 

tube-side: 
section i: 


= 0 

de 

sections j=2 to J=H: 


dtj 

de 


* A;J-1 ) - Vi (t j-t )/DZ 


A. 

fij di Cpj 


(4-1-7) 


tube- wall : 

sections J-1 to j=Jl: 


* 'Bi (Tj-Twj) - Cj (Twj-tj) 


de 


A hoi do 4 hi.i di 

® J * fwCpwTdo2Idl27, ^ ” fwCpw(doii-dii^; 


(4-1-8) 


shell -side: 

sections J=1 to J=n-1: 


= -voj (^^5i=^) - Dj (Tj-Twj) 


DZ 


TN, _ 4 ho.i do 

^ ^ojCpoj (ds2-do2 ) 


(4-1-9) 


where ds: shell inside diameter 
section M: 




1 
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The resulting equations are then solved by one of the 
numerical-integration methods discussed in Chapter 3. 


q.1.2 Example and Results 


A computer program was written to simulate the dynamic 
character! Stic of the heat exchanger. The program, which uses the Euler 
method with 0.018 sec. time intervals, was run on an IEM/370 to 





•*y ' 





I 


simulate the response of a step change (50 ®K) in the shell-side input 
temperature. Table 4-1 shows the input data and the results are given 
in Figure 4-3 and 4-4. The time scale in Figure 4-3 is 100 times of 
that in 4-4, which shows a more detailed change at the shell-side. 

As expected, the rate of the steadily increasing temperature at 
tube-side was much slower than that at the tube-wall and shell-side. 



The final outlet temperature at the shell-side increased by 33 ®K and was 
66% of the inlet step change. Finally, although this is a nonlinear 
system, the response of the shell-side is similar to a first-order 
response and the tube-side and tube-wall are similar to second-order 
response with a damping factor greater than 1. 


Definition 



Dimension Initial 
value 

Unit 

3ss>: 

4 

0.63 

g-mole/s 


6 

49.1 

g-mole/s 


7 

130.7 

g-mole/s 

n 


1.36 

atm 

TCO 


300 

K 

ss 


1.83 

m 

31 


0.045 

m 

32 


O.05O8 

m 

33 


0.0732 

m 

n 


55 


r. 


6 


an* 


8027 . 17 

kg/m 

CPV 


0.502416 kJ/kg-K 



20.76882 

j/m-s-K 

Xtso 

1 

1.49 

g~mole/s 


2 

1.37 

g^mole/s 


3 

16.92 

g*mole/s 


4 

14.28 

g-mole/s 


5 

70.77 

g-mole/ s 


6 

0.96 

g-mole/s 

P3 


4.26 

atm 

?HZ 


531 

K 


Flow rate of H2O in tubeside 
Plow rate of N2 in tubeside 
Flow rate of O2 in tubeside 
Pressure of tubeside 
Temperature of tubeside inlet 
Length of heat exchanger 
Inside diameter of tube 
Outside diameter of tube 
Inside diauneter of shell 
Number of tubes 

Number of finite-difference sections 
Density of wall 
Heat capacity of wall 
Thermal conductivity of wall 
Flow rate of CH^ in shellside 
Plow rate of CO in shellside 
Flow rate of CO2 in shellside 
Flow rate of H2O in shellside 
Plow rate of H2 in shellside 
Flow rate of N2 in shellside 
Pressure of shellside 
Temperature of shellside inlet 



Table 4-1 Input data of dynamic simulation of heat exchanger 
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4.2 Shift Converters 


The shift converters discussed in Section 2.1.2 is remodeled here 
for dynamic analysis. The reactor, which was treated as a fix-bed 
reactor, has no provison for removing the heat of reaction, such as a 
Jacket, so it operates with an adiabatic process. The heat generated by 
the reaction, therefore, heats up the gas as it flows through the 
reactor. The gas then exists at a higher temperature than the inlet 
feed temperature. ^1 the assumptions described in Section 2.1.2 are 
also applicable here. 


4.2.1 Ilathematical tlodel 


The dynamic modeling of the shift converters must take into account 
the energy balances of reacting gas and catalyst, and the material 
balances of reacting gas. 


Similar to the distributed model at steady state, the mass balance 


3 c/36 = -V 3 C/3Z - ra*5g/£ 

and the terms are the same as those in the steady state model. 


(4-2-1) 




The two energy balazice equations are 

reacting fluid: Cp p ^t/ae = -Cppv bt/az + <-Ah^) (ra» /e> 

♦ hcAc (Tp“t> 
catalyst: Cpc Z Tp/dO s'hckc ft-Tp> 

'diere he: particle-fluid beat transfer coefficient, 

Ac: specific .surface area of packing, n*/m* 


(4-2-2) 

(4-2-3) 
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Cpc: heat capacity of catalyst, J/kg“K 
Tp: temperature of catalyst,K 

The Initial conditions are the distributed model's steady state 
solutions, where the catalyst temperature is assumed to be the same as 
the fluid temperature, and the boundary conditions are 


dt inlet 

da 




dc inlet 
dS 


(4-2-5) 


The only parameter which has not been discussed is he, the 
particle-fluid heat transfer coefficient. From Ref. 32 he can be 


calculated by 


t Jh s 0.395/Re 
djh = 0. 725/Re 


30 ^ Re ^ 10 
2 ^ Re ^ 30 


where Jht Colburn analogy factor of heat transfer, and 
Jh = h /V p Cp 


(4-2-6) 


(4-2-7) 


(4-2-8) 


h: heat transfer coefficient, J/m*-s-K 


v: velocity of fluid, m/s 

p: molar density of fluid, g-mole/m* 

Cp: molar heat capacity, J/g-mole-K 


Pr: Prandtl number 


Re: Reynold number 
: £ : void fraction 


P^fferential-Dlfference Equations 



The distance varying terms were changed to finite-difference form 
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(back-difference). The final set of differential-difference equations 
for the reacting fluid are: 

Section 1: 


dc, 

de 


= 0 


dt 

de 


= 0 


(^-2-9) 


sections j=2 to j=h: 


dcj 

de 


—J-l 


( 4 - 2 - 10 ) 




dt 

de 


"^3-1 ^ BZ ^ * Cp j-1 p 0-1 




Cpj pj' 


Ac(Tpj-tj) 


( 4 - 2 - 11 ) 


cavalyst : 

sections j=l to j=K: : 


m 




d^Pj hcJ ftc , 

dS = Cpc-e>Q - Tpj ; (q-Z-lZ) 

These simultaneous differential -difference equations^ (4-2-9) through 

(4-2-12), were solved to obtain the temperature and concentrations of 

the reacting fluid under transient conditions. 


^■2.2 Exanole and Results 


The computer program, which utilizes the Euler method with a 
0.00072 sec. time interval, was run on an IEIt/370 computer to simulate 
the transient responses of the high temperature shift converter. Table 
4-2 shows the example's input data and the results are shown in Figure 


4-5. 




Figure 4-5 shows ^hat the outlet temperature did not change but the 






mm 





Dimension 

Initial 

value 

Unit 


1 

1.52 

g^mole/s 

2 

10.93 

g-mole/s 


3 

7.53 

g-mole/ s 


4 

23.95 

g-mole/ s 


5 

62.91 

g-mole/s 


6 

1.05 

g-mole/s 

n 


1.5 

atm 



696.2 

K 

s 


1.8288 

m 

y 


0.09144 

m 

r 


92 


f 


6 


UJ 2 X 


964.57 

m2/m3 



0.879 

kj/kg-K 


1281.5 


0.469 

0.00505 


i 


Definition 


kg catalyst/m5 
bed 


Inlet flow rate of CE4 
Inlet flow rate of CO 
Inlet flow rate of CO2 
Inlet flow rate of H2O 
Inlet flow rate of H2 
Inlet flow rate of N2 
Inlet pressure 
Inlet temperature 

Height of high temperature shift converter 
Diameter of high temperature shift 
converter 
Number of tubes 

Number of finite difference sections 
Specific surface area of packing 
Heat capacity of catalyst 
Density of packing 


Void fraction 
Diameter of catalyst 


Table 4-2 Input data of dynamic simulation of high temperature 
shift converter 
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1^2 


:;onversion rate Ss increased. The reason why this situation occurs can 
be observed from Equation 4-2-2, v^iere the time varyir^ temperature 
depends on three terms which are the material . transport, the reaction, 
and the heat transfer terms. The reaction term is directly proportional 
to temperature ( exothermic reaction) while the other two terms are 
Inversely proportional. The temperature increase at the inlet is 
balanced inside the reactor. Actually, most of the increased heat 
energy is used to increase the temperature of the catalyst in this 
example. 

Because the rate expression in shift converter is an expotential 
function of temperature, the accumulation term is very sensitive to the 
temperature and the stable time interval is small. The small time 
interval results in a larger trunction error and a longer running time 
for the computation. A method for improving this situation is discussed 
at the end of this chapter. 
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4.3 Reformer 


The combustion gas heated steam reformer vas discussed in Section 

2.1.3 where it was modeled as a heat exchanging fixed-bed reactor. 
Equations were derived to determine the temperature and concentration 
gradient as a function of location in the reformer only. In order to 
eliminate time as an additional independent variable, steady-state 
operation was assumed. The more general case of unsteady-state 
operation is discussed here, with the the dependent temperature and 
concentration variables as functions of the independent variables, time 
and distance. The reformer is shown in Figure 2-E. Because there is 
more than one Independent variable, the relationship between 
temperature, flow rate of the components, time, and distance can be 
stated as a partial differential equation. However, since the objective 
is to reduce the mathematical model to a set of simultaneous 
differential-difference equations (i.e., only one independent variable, 
time), the system will be "finite differenced" with respect to the 
location. 

^.3.1 Hathematical Kodel 

A detailed mathematical model for studying the dynamic response of 
the reformer was developed. All of the assumptions stated in Section 

2 . 1.3 are applicable here. 





c = P Ye 
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Z> Z 


_ -d(T^e) 4 
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where F = molar flow rate inside the reactor, g-mole/hr 
di = inside diameter, m 
Equation (4-3-3) is replaced by 


(4-5-4) 


■a?Ye 

02 


ra'ifB 


4 


0 

ae 


(e^e)^ 


(4-3-5) 


The conversion, xl, can be replaced by the ir.ole-f taction of methane: 


Ye 


Fo+2xi Pi 


VI - Fi~Fo Ye 

" P142Fi Ye (4-5-6) 


where Fo = initial total molar flow rate inside reactor, g-mole/hr 
F^ = initial methane flow rate inside reactor, g-mole/hr 
The equations for the molar flow rate and density as a function of 

9 

conversion are given by: 

F = Fo + 2x1 Pi (4-5-7) 



where P = pressure, atm 
R = gas constant 

t;* tube-side , tempera'^^e, K . . ■ _ : 


(4-5-8) 


Ejiergy Ealanse 


lAiS 


Four onorey balances are required tor the reformer: one for the 
reformer gases (tube-side), one for the tube-wall, one for the 
combustion gases (shell-side), and one tor the catalyst. The reformer 
gas balance includes its own sensible heat change, reaction enthalpies, 
heat transfer from the hotter tube-wall, heat transfer to the oatalyst 
particles, and the accumulation term, which are as follows: 


a FCpt 

a z 


P (-di H) H + hi A di(Tw-t) + hcAc 


dZ 




ae 


( 4 - 3 - 9 ) 


where Cp = heat capacity on the tube-side, J/g-mole-K 
X = tube-side ter.perature, K 

he a exterr.al surface area of particle per volume of catalyst 
bed, ir.VjT.’ 

= tube-wall temperature, K 
Tp = catalyst temperature, K 

hi s inside heat- transfer coefficient, J/hr-m-K 

he s particle-fluid heat transfer coefficient, O/hr-tr-Y. 

AH = heat of reaction, J/g-mole 

In the heat of reaction term, the component molar flow rate, F, is for 
that component corresponding to the conversion x. Therei-ore, fo* 
reactions at hand . . • 1 ; 

ZP(-d^H)||= PiC-AHi)|f + (P5«lPl) 

o ^ (4-3-10) 


where JiH. = the demethanat ion reaction errthalpy, J/g-mole CH4 


1^7 


= the water shift enthialpy, J/g-mole CO 
The enthalpies, , are evaluated at the average reformer temperature. 
Negative values indicate an exothermic reaction. 


Shell-side and tube-wall equations for the energy balance are as 
follows: 


pc Ao Cpo = -^0 VO Ao Cpo - ho^do(T-Tw) 


(^w(do2-di2) Cpw 1^ = 4 ho do(T-Tw) - 4 hi di(Tw-t) 


where = Fluid density ?t shell-side, g-nole/m* 

Ao = shell-side flow area, m* 

Cpo = heat capacity on shell-side, J/g-mole-K 
T = temperature of shell-side, K 
Vo = shell-side velocity, m/hr 

ho = outside heat-transfer coefficient, J/hr-m*-K 
Tw = tube-wall temperature, K 
do = outside diameter, m 
pw = density of tube-wall, kg/ra’ 

Cpw = heat capacity of tube-wall, J/kg-K 


The energy balance equation for the catalyst is 


-eg = he Ac(t-!I^) 


'Acre Cpc = heat capacity of caralyst; J/kg-K 


The initial conditions for these equations were taken as the steady 
state solutions and the boundary conditions were the new inlet 
variables. 

Solving the Kathematical Kodels 

Expandin'; Eq. (<l-3-5) and Eq. (41-3-9) to obtain; 
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ra* <£3 
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(4-3-14) 





= - *2; r(--iH) || + hiXdi(TM-t) 


♦ he Ac 


4 


(Tp-t) 


(4-3-15) 


Expanding oP/3 0 by using Eq. (41-3-6) and (4J-3-8) and the following 
identity: 


af _ f /gxi^ 

^ 0 “ V / 'aYo4 V fl 


'‘^Ye' '•a 6 


) 


making the necessary substitutions into this identity, we will obtain: 


^ sr 2P Fi 

(4-3-16) 


|2l _ -(Pvf2Fi Ye)'Po - (Fi~Ye Fo)-2-Pi - 

0 *® (P 1 + 2 P 1 ye')ii ( 4 - 3 - 17 ) 


. g_P Pi / -(PH-2P1 Ye)-Po - (Pi-Ye Po)-2-Pi n 9 Ye 
fi t Po ' (p^+2Pi Ye;^ ' a e 


(4-3-18) 


1<59 


In order to find solutions for the variables in the example for a 
tubular reactor, numerical data for all the variables at time 9 is 
needed. In using the finite difference technique, the reactor is 
considered to have H sections (see Figure 2-7 ) and values for the 

variables at (0 + A 0) are obtained. The final set of 

differential-difference equations are as follows: 
tube-side: 

sections J=1 to j=M, define 

2 P.1 Fi / -(F1^^2F1 Ye.iVFo - ( Fi~Ye.i FoV2-7ia 
' R tj Fo ^ (Fi+2Pi Yej)2 

Cj = D ( Pd+Yej Aj) 

Rtj Fo ^ 

Pj = Fo + Fi 2 xij 

Gj = D Cpo 


D = 


^dl2^ 


(4 


section l: 


liL = 

de 


0 , 


d Yel 

d e 


= 0 , 


d Pi 

d e 


= 0 


0 


•ections j=2 to 

= -(F.1 Ye.1 - Fi -1 Yei- 1 ) _ ra'i-iesD 
**0 DZ Cj Co € 


(4- 


iLl - ,./d Yejs 

« 0. r- 


(4 
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define 
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.. (dP.:s 


D Cpj 
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dtj 

de 


1 r / CT3.1 ?,1 t.i - CT3-i-1 F.i-1 ^ 

= 1“ ^ DZ 

(?1 (-AE 1 J. 1 ) * (j'5 + =‘10-1 *'i)(-‘‘'^20-i) ■" 

hcO Ac |<TpO-tO) + .--iO-l(Adi) (T«3-1-tj-l) ' ^o] (4-;-2i) 


shell-side: 
sections J=1 to 


d?.i -vo-i (T-i- 1.1 - Ti) „ /„ 

d6 = 

where 


.. . _ ^ hoi do 
"0 “ (-^oj Ao Cpoj 

section H: 


( 4 - 5 - 25 ) 


ir = ° (4-5-26) 

tube-wall: 

sections j=l to j=tl: 


de^ = Oj (Tj - Twj) - Qj (Twj - tj) 

- 4 hOi do . ' - . • c- . " . C--. 

^ Pw Cpw (do2-di2) 

t ^ '= ~ 4 hi-i di ' ' " " ' - ' - - 

. Q w Cpw (do2-di2 ) , . -■ ( 4 - 3 - 27 ) 

ORiGIMAL P^&IS 
PODI^-t^UAUtY 


cstalyst; 
%ec^ions j=l 


(4-3-2 



he-? Ac 
Cpc €3 


(tj - Tpj) 


il. 3.2 Example and Results 


The corresponding computer program, which use the Euler method with 
a 0.0036 sec. time interval, was run on IBK/370 computer to simulate the 
dynamic response of the reformer. The initial disturbance was a 100 
step increase in the inlet combustion gas temperature. Tables 4-3 and 
4-4 present the input data and the output of the transient response, 
respectively. 


Note that the outlet temperature on the combustion gas side 
increased by 16 *K while inlet temperature increase was 100 *K , and the 
outlet temperature of the reforming gas did not change much. Actually, 
the reforming gas outlet temperature would gradually increase with time. 
The increased energy is partly going into heating the wall (then 
reforming gas and catalyst) and partly simply "wasted” with stream flow. 
Here the competition between heat transfer and flow rate can been seen. 
A more complex condition in the reforming gas, occurs when the reaction 
rate is the third competitor (for an endothermal reaction) and more than 
one direction for heat transfer Is allowed. 


Again, the small value for the time interval is a problem, which 
required much computation -time. In the following chapter, methods for 
Improving this condition are discussed. 


Befini-tion 


••triable Biriensior: Initial Unit 


value 


0»0375 g-nole/s Input CH^ flow rate 
3 Stear. to carbon ratio 


FX 


1.5 

atm 

Inlet Cr*4 pressure 

t:c 


811.1 

K 

Inlet temperature 

ZF. 


1.524 

n. 

Height of reformer 



0.103 

m 

Outside diameter of regenerative tube 

D2 


0.1315 

rr. 

Inside diameter of reforming tube 

B 3 


0.1414 

m 

Outside diameter of reforming tube 

KC 


10400 


Rate constant of demethanation reaction 

£A 


83736 

j/mol e 

Activity energy of denjethanation 





reaction 

HHC '3 


1281.48 

kg/m3 

Density of packing in reformer 

LFS 


0.469 


Void fraction in reformer 

S 


0.229 

m 

V.'idth of combustion gas square duct 

Di' 


0.005 

m 

Diameter of catalyst 

K 


23 


Number of finite difference sections 

Rr>’ 


8027 . 17 

kg/m3 

Density of vail 

CFW 


0.6113 

kJ/kg-K 

Heat capacity of wall 

SASIA 


669.29 

m 2 /ni 3 

Specific surface area of catalyst 

crc 


1.026 

kJ/kg-K 

Heat capacity of catalyst 

LK 3 G 

5 

0.042 

g-mole/s 

Inlet C02 flow rate in combustion gas 


4 

0.036 

g-mole/s 

Inlet H 2 O flow rate in combustion gas 


6 

0.109 

g~mcle/s 

Inlet 3v2 flow rate in combustion gas 

F 3 


1.5 

atm 

Inlet combustion gas pres s\are 

TH 2 


1811,1 

K 

Inlet combustion gas temperature 


Table 4-3 Input data for dynamic simulation of reformer 
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-g. 


■hr. 


0.00000 
0.50000D-04 
O.lOOOOD-03 
0.15000D-03 
0.20000D-03 
0.25000D-03 
0.30000D-03 
0.35000D-03 
0.40000D-03 
0.45000D-03 
0. 500000-03 
0. 550000-03 
0.600000-03 
0.650000-03 
0.700000-03 
0.750000-03 
0.800000-03 
0.850000-03 
0.900000-03 
0.950000-03 
0.100000-02 
0.105000-02 
0.110000-02 
0.115000-02 
0.120000-02 
0.125000-02 
0.130000-02 
0.135000-02 
0.140000-02 
0.145000-02 
0.150000-02 
0.155000-02 
0.160000-02 
0.165000-02 
0.170000-02 
0.175000-02 
0.180000-02 
0» 185000-02 
0.190000-02 
0.195000-02 
0.200000-02 


Outlet temperature Outlet Last section Last section 

of reforming gas temperature wall catalyst 

K of combustion temperature temperature 




gas, K 


K 


K 


0.951950 

03 

0.10348D 

04 

0.99917D 

03 

0.95195D 

03 

0.951950 

03 

0.10348D 

04 

0.99919D 

03 

0.95195D 

03 

0.951950 

03 

0.10348D 

04 

0.99922H 

03 

0.95195D 

03 

0.951950 

03 

0.10348D 

04 

0.99924D 

03 

0.9519511 

03 

0.951950 

03 

0.10348D 

04 

0.99926D 

03 

0.95195D 

03 

0.951960 

03 

0.10348D 

04 

0.99928D 

03 

0.95195D 

03 

0.951960 

03 

0.10348D 

04 

0.99930D 

03 

0.95195D 

03 

0.951960 

03 

0.10348D 

04 

0.99932D 

03 

0.95195D 

03 

0.951960 

03 

0.10348D 

04 

0.99934D 

03 

0.95195D 

03 

0.951970 

03 

0.10348D 

04 

0.99937D 

03 

0.95195D 

03 

0.951970 

03 

0.10348D 

04 

0.99939D 

03 

0.95195D 

03 

0.951970 

03 

0.10348D 

04 

0.99941D 

03 

0.95195D 

03 

0.951970 

03 

0.10349II 

04 

0.99943D 

03 

0.95 19511 

03 

0.951980 

03 

0.10350D 

04 

0.99945D 

03 

0.95195D 

03 

0.951980 

03 

0.10353D 

04 

0.99947D 

03 

0.95195D 

03 

0.951980 

03 

0.10358D 

04 

0.99949D 

03 

0.95195D 

03 

0.951990 

03 

0.10366D 

04 

0.99951D 

03 

0.95195D 

03 

0.951990 

03 

0.10375D 

04 

0.99953D 

03 

0.95196D 

03 

0.951990 

03 

0.10388D 

04 

0.99956D 

03 

0.95196H 

03 

0.952000 

03 

0.10402D 

04 

0.99958D 

03 

0.95196D 

03 

0.952000 

03 

0.10417D 

04 

0.99960D 

03 

0.95196D 

03 

0.952000 

03 

0.10433D 

04 

0.99962D 

03 

0.95 19611 

03 

0.952010 

03 

0.10448D 

04 

0.99964D 

03 

0.95196D 

03 

0.952010 

03 

0.10462D 

04 

0.99966D 

03 

0.95196D 

03 

0.952010 

03 

0.10474D 

04 

0.99968D 

03 

0.95197D 

03 

0.952020 

03 

0.10485D 

04 

0.99970D 

03 

0.95197D 

03 

0.952020 

03 

0.10494D 

04 

0.99972D 

03 

0.95197D 

03 

0.952030 

03 

0.10501D 

04 

0.99974D 

03 

0.95197D 

03 

0.952030 

03 

0.10506D 

04 

0.99976D 

03 

0.95197D 

03 

0.952030 

03 

0.10510D 

04 

0.99978D 

03 

0.95198D 

03 

0.952040 

03 

0.10513D 

04 

0.999BOD 

03 

0.95198D 

03 

0.952040 

03 

0.10516D 

04 

0.99982D 

03 

0.95198D 

03 

0.952040 

03 

0.10517D 

04 

0.99984D 

03 

0.95198D 

03 

0.952050 

03 

0.10518D 

04 

0.99986D 

03 

0.95198D 

03 

0.952050 

03 

0.10519D 

04 

0.9998BD 

03 

0.95199D 

03 

0.952060 

03 

0.10519D 

04 

0.99990D 

03 

0.95199D 

03 

0.952060 

03 

0.105200 

04 

0.99992D 

03 

0.95199D 

03 

0.952070 

03 

0.105200 

04 

0.99994D 

03 

0.95199D 

03 

0.952070 

03 

0.105200 

04 

0.99996D 

03 

0.95200D 03 

0.952070 

03 

0.105200 

04 

0.99998D 

03 

0.95200D 

03 

0.952080 

03 

0.105200 

04 

O.'IOOOOD 

04 

0.95200D 

03 


Table 4-4 Output of the dynamic simulation of the reformer 
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Discussion 


The objective of this chapter was to exanine and develop a computer 


sinuilation for the dynamic characteristics of components in the fuel 


process subsystem. The main program and subroutines for each componeni 


will be modified to be a group of subroutines in the final simulation of 


whole system, and some of the duplicated subroutines will be omitted. 


The examples illustrated in this chapter were not concerned with 


the flow rates which were assumed to be constant. However this does 


affect the transient responses due to transpotation lag in the pipe. 


This factor is considered in the next chapter. 


The selection of an appropriate time interval in solving difference 


equations is a problem in digital simulation. Large values for time 


Intervals will "blow-out" the computation. As the calculated value 


approaches infinity. However, smaller time interval values require more 


computation time to iterate, and increasing truncation error occurs in 


the computation. One possible solution to this problem may be the 
Runge-Kutta-Kerson integration method (Section 3-<l), which automatically 
•elects the proper time interval. However, this method does not work 
**•11 with a reaction term that is expontential with temperature change. 
^ alternate method, which is discussed in the following paragraph, 
places the attention on the mathematical model. 



For the heat exchanger, it can be observed from Equation 4J-1-8 that 

'lie D2/V 


I 


term (the time needed to travel on section) determines the time 



1B5 


Interval (H). The change in temperature of the previous section will be 
greater if H>DZ/v. Therefore, DZ/v is the maximum H value, especially 
when the transport rate is larger than heat transfer rate, which is 
usually the case. It is not the same condition in the catalytic 
reactor, where the maximum H is determined by the reaction rate which 
is much faster than the other two. If we add one more assumption to the 
simulation, such that no temperature difference exists between the 
reacting gas and the catalyst, then Cp in the accumulation term (see 
Equation <1-3-9) is the composite heat capacity of gas and catalyst, 
which is much larger than that of the gas only. Subsequently, this 
assumption will decrease the amplification, and thus increase the stable 
time interval. In the example for transient simulation of the reformer, 
this additional assumption increases the time interval from 0.0036 sec. 
to 0.05 sec. The simulation of the whole system discussed in Chapter 6, 
includes this assumption. 


If 


CHAPTER 5 


DYNAMIC MODELING OF THE FUEL CELL STACK SUBSYSTEM 




As stated in Chapter 2, within the fuel cell stack, hydrogen and 
oxygen react to continuously produce DC electricity, waste heat, and 
steam. 


Generally, the operation of the PAFC system during a load change is 
to control the voltage level to obtain the required power level, A 
change in the voltage level will affect the current density on the cell 
plates. The relationship between the current density, voltage, and 
power is obtained from Ohm's law and the current-voltage characteristics 
stated in Chapter 2. One typical result is shown in Figure 2-27 for 
fixed operating conditions. 








KM 


Hhen the current density is changed, the inlet flow rates should 
also be changed to meet the requirements. Therefore, the variables 
that can be manipulated for control of the fuel cell stack subsystem 
are the input flow rates of fuel (hydrogen) and oxident (oxygen or air). 


In addition to the assumptions made for calculating the steady 
*^te conditions in the stack, there are several other assumptions made 
for system simulation during the load change period: 

f • The temperature and composition of both of the reservoirs (fuel and 
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air) outlets (see Figure 5-1) do not change. However, the flow rates 
are changed in proportion to the DC change. Also, the inlet temperature 
of the process air does not change. 

2. The pressure in both reservoirs is slightly above the maximum (full 
load) operating pressure. 

3. When changing the fuel cell operating temperature and pressure, a 
line of constant phosphoric acid concentration is to be followed. 
Figure 5-2 shows the total pressure of phosphoric acid solutions as a 
function of temperature. This criteria produces a significant increase 
in the estimated cell life at 25% of the rated power (Ref. 27 ). 

4. The utilization ratios of both fuel and process air are the same. In 
addition, the stoichiometric factor for the cooling air is hold 
constant . 

5. Because the reaction rate (including the reacting gas diffusion and 
ion migration) is faster than the heat transfer rate and the component 
flow rate, the accumulation of components resulting from the reaction is 
neglected. 



Page Intentionally Left Blank 




Equilibrium HjPO^ concentration 
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5.1 Design of The Subsystem for Load Change 


A preliminary process flow diagram, for control of the load change, 
is shown in Figure 5-1 . The basic control steps for a load change in 
this subsystem are as follows: 

1. The flow control valves (FCV-1 and FCV-2) regulate the flow of fuel 
and process air in direct proportion to the DC current. 

2. The cooling air flow control valve (FCV-3) regulates the cooling air 
flow rate in proportion to the flow rate of the process air. 

3. The governor control valve (RPM) regulates the speed of the air 
compressor (air flow rate is proportional to the speed). At steady 
state conditions the inlet air flow rate is the same as the process air 
flow rate required in the stack; however, the air flow rate depends on 
the pressure level of the reservoir in transient conditions. 

4. The pressure levels of the fuel and process air are regulated by the 
pressure control valves (PCV-l and PCV-2 respectively), while the 
cooling air pressure level is the same as that of the air reservoir and 
is complemented by the circulator through speed or pitch control. 


The advantages of placing the two reservoirs before the fuel cell 
stack are as follows: 

1. The flow rates can be rapidly changed because the independence of 
fuel processing subsystem. 

2. The fuel cell stack subsystem can operate for a period of time when 
the air compressor is not working. 
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5.2 Hathematical Wodel 




The mathematical model for evaluating the transient temperature 
distribution in the stack is similar to the model used to simulate the 
steady state temperature distribution. {Additional accumulation terms 
are used in the transient model to represent the unsteady state 
situation. {According to assumption (5), there is no accumulation of 
components caused by the reaction. The model for the current density 
distribution, which interprets the mass balance and was discussed in 
Section 2.2.2, is applicable for calculation under transient conditions, 
however, the inlet flow rates and the temperature at each grid section 
arc functions of time. 


Considering one strip of the PIAFC stack (see Figure 2-17), we can 
write the energy balance for the plates in the form of 
differential-difference equations. 

(1) Engergy balance In the stack 
(a) cooling plate 
section 1: 




«=01Cc(i 

Pc 


rc(i Vtco; 
A y 






sections 2 to M-1; 




T(2,3) - (Tc(j)-Tc(o-I)) =pc cc t* ^ 
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pc aY 




= ^c cc t' 


dT(l.M) 

de 


(5-2-5) 


(b> cell plate No. 1 


section 1 * 


— ■7’ AV*-!*-/ ' 

+ ^ '^^5,1) 

AX2 


mpCia: 

Pp 




T(1,1) 

+ [v*(l)-v(l)] 1 ( 1 , 1 ) 


= fp cp 


A <iT(2.ll 

^ de 


(5-2-4) 


sections 2 to M-lS 


Kyt ( 2Kyt ^ K2- + 1^) T(2,o) +-§Z T(2,j-l) + ■ 

ISOj. T(2,j+1) - (-^+AXi ax2 

"L SEil^j^Cpi:^ (Tp(1,3)-T^^^ 

+ HjT(5,j) - -pj^y 


section H: 


- (^ S * St ^ 

im(2.Mlgo(2.Wl (Tp(l,M)-Tp(l.»-J^) 

“ PpA.y r- ' ,vrr- ■■■ C-':^ ■• 

4Tl2a»l; , 


(5-2-6) 




I 


I 


(c) cell plate No. 2 


section 1: 




mp(3.l)CT 
Pp A y 


(Tp(2,l)-Tpo) +[(V«(2)-V(2)) 1(2,1) 


= Pp=P * ^ 


sections 2 to W-1: 


^T(3,j+l) - T(3,J) +^T(3,j-l) +^T(2,i; 


- s 

.p,cp t 


+[(V*(2)-V(2))j 1(2,3) 


section M: 


^ T(5,,^1) - (^ . ^) T(3,M) T(2,«) . g. T(4.„) 

■ (Tp(2,M)-Tp(2,M-l)) + j(V»(2)-V(2))j I(2,K) 


“PpPPt^ 


plate No. 3 (symmetry plate) 

*«ctlon 1: 

. . ; i T - ■ 

^ T(4,2) - + 

■ (i^0.i‘)-ii>o)*4|(yi(5pv(3y|i^^ 

=e, =P t 


( 5 - 2 - 10 ) 


16^ 


sections 2 to M-1: 




gpn.ji 

Pp^y 

Pp cp t 


(Tp(5,j)-Tp(3,j-1)) +[(v*(3)-v(3))]l(3,j) 


section H: 


T(4.M-1) T(4.M) * Ig T{5,«) 

_ rg)(3tM) Cp(?,Ml (Tp(3,M)-Tp(3,M-l)) +f(V*(3)-V(3) )] l(3, 0 ) 
Pp A y *• 
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Initial Conditions : T(x,y,0)=Ti (x,y) — steady state distribution 

where 

Cc(j) : molar heat capacity of cooling air in element (i/J)» J/g-mole-K 
Cp(i, j) : molar heat capacity of process air in element (i»j), 

J/g-mole-K 

I(i,j) : current density of element <i,j>. A/cm* 

kx : thermal conductivity along stack direction, J/hr-m-K 

ky : thermal conductivity in flow direction, J/hr-m-K 

M : number of elements along air flow direction 

mc(j) : molar flow rate of cooling air in element (i,J), 

g-mole/hr-channel 

■ip(i, J) : molar flow rate of process air in element (i, J>, 

, - ■ ' c;;; ' ■; g”mole/hr-channel 

Pc : pitch of cooling cluinr^l,^,m ,, j. . 

Pp : pitch of process air channel, m 

t : thickness of cell plate including matrix, bipolar plate. 


and electrodes, m 


t' : thickness of coolins plate, m 
T(i,j> : temperature of cell plate, K 
i = 1 cooling plate 
i = 2 cell plate No. 1 
i = 3 cell plate No. 2 
i = 4 cell plate No. 3 

j : sections defined along air flow direction 
TcO : inlet temperature of cooling air, K 
Tc(j) : temperature of cooling air, K 
TpO : inlet temperature of process air, K 
Tp(i,j) : temperature of process air, K 
i s 1 plate No. 1 
1=2 plate No. 2 
i = 3 plate No. 3 

V<i) : operating voltage of cell plate i, V 

V*(i) : ideal voltage of plate 1 under operating conditions, V 
xl : t/2 + t'/2, m 
x2 : t, m 

y : finite difference distance along air flow direction, m 
^p : density of cell plate, kg/m^ 

^c : density of cooling plate, kg/m* 
e : time, hr 

cc : specific heat, of cooling plate, "Vr^,T--Sor- ' 

CTp : specific heat, of cel^ J/kg-K„ th- 
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(2) Convective heat transfer (energy balance) in the channels (both 
process and cooling air) 

^ RoC iJ/de = hs(Tw-T)-mC 3T/Sx 
B.C. X = 0 T = Tj-l (0) 

I.C. 0 = 0 T = T' 

where 


Tw 

Tj_i 

■ti 

h : heat transfer coefficient, J/m*-hr-K 
s : perimeter of the channel, m 
m : molar flow rate per channel, g-mole/hr 
c : molar heat capacity, J/g-mole-K 
Tw : temperature of cell plate, K 
T : temperature of air, K 
X : distance from the edge, m 
9 I molar density, g-mole/m* 

Ao : cross area per channel, m‘ 

6 : time, hr 

' : steady state solutions 

The equations were solved by the Laplace Transform method, and the 
solutions were changed to difference form with the properly chosen time 
intervals (detailed discussion is in Appendix 4). 




for process 


air 


, K TjCH) = To-i(0) exp (^) + Twj ( 1-exp C^")) 

=2:-: Tj(2H)= Tj_i(H) exp (=|^) + Tw j( 1-exp (“ y—) ) 

etc. 


for cooling air 


0 = E Tj(H) = Tj-l(H) exp + Tw^Cl-expC-^ — )) 

0 =2H Tj(2H)= Tj_i(2H)exp + Twj(l-exp(^|^) ) 


etc. 

where 


A = 


hs 

PAC 


B = 


m 

^A 


Twj = (Tw(e) + Tw(6+H))/2 
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5.3 Computer Program 

The procedure for the simulation of the transient respones of the 
current density and temperature distributions on the cell plates, when 
the system load decreases, are as follows: 

1 . Obtain the steady state solutions. 

2. From time=0, decrease the load at a fixed ramping rate until the 
required load ratio (the minumum ratio is 25%) is reached. 

3. Calculate the corresponding operating pressure to keep the phosphoric 
acid (electrolyte) concentration constant (see Figure 5-2). 

4. Estimate the mean current density and operating voltage at the new 
operating conditions to fulfill the new load. Two typical figures 
(Figure 2-27 and Figure 5-4 ) show the relationships between the power, 
current density, and operating voltage. 

5. Calculate the flow rates of the fuel, the process air, and the 
cooling air at each position from the new mean current density. 

6. Estimate the current density distribution from the new flow rates and 
pressures. 

7. Calculate the new temperature distribution on the cell plate. 

8. Calculate the new temperatures of the gases at each position. 

9. Stop the procedure when the final (new) steady state is reached. 
Otherwise, increase the time by a time interval and go back to step 2. 



This procedure was developed into a Fortran computer program. The 
final set of differential-difference equations described in previous 
section was solved by ' one of the integration methods discussed in 
Section 3.4. In Figure 5-3, a simplified flowchart of the computer 




IFUC=0 
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program, where the status of various variables (current density, flow 
rate, temperature of plate and gases, pressure, voltage, and power 
output) at different times (0 and0+A0) is shown. Figure 5-3 also 
expresses the function of each calling subroutine. 
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5.^ Results and Discussion 

The program, which uses the Euler method (Section 3.4) with a time 
interval of 0.108 seconds, was run on an IBlt 370 computer and the 
results are shown in Figures 5-5 through 5-16. Using a 12.5% per second 
ramping rate, the load decrease (from full to 50%) can almost be treated 
like a step change. Table 5-1 lists the input data in SI units. The 
cooling configuration is a "straight" type (see Section 2. 2. 3. 3) and the 
inlet flow temperature is close to the average operating temperature. 

Figures 5-5 through 5-10 show the numeric distributions of the 
temperature and current density in each finite-difference cell on the 
symmetric plate (plate no. 3) at the initial steady state, the transient 
state (10.152 sec.), and the final steady state (280 sec.). The plots 
of the temperature distribution for the three cases are shown in Figure 
5-11 and 5-12. Finally the mean values of the temperature, voltage, 
and current density at special points were plotted vs time, with a small 
time scale and a large time scale. The results are shown in Figures 
5-13 through 5-16. From the results of the small time scale plot, the 
changes in the early transient state, which decide the stable ramping 
rate and dead time for control, can be seen. The plots with the larger 
scale provide Information concerning the transient time period, the 
system response between the two steady state periods, and the final 
steady state conditions. 

As expected, the temperature responds slower than the current 
density which is dependent upon the inlet flow rates. The damping 
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Figure c- 10 Final steady state current density distribution (200 sec.) 
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Definition 


Variable 

name 

Dimension Initial 
value 

Unit 

XN 

0.41 

m 

YN 

0.28 

m 

XDNSO 

3250 

A/cm^ 

UTA 

0.5 


UTH 

0.75 


POPC 

3.4 

atm 

POP 

3.4 

atm 

TKA 

443 

K 

WFD 

0.001016 

m 

WF\/ 

0.003048 

m 

NCC 

30 


'WE 

0.001016 

m 

TKF 

450 

K 

4. 

0.003302 

m 

NX 

5 


WAD 

0.001016 

m 

WAW 

0.003048 

m 

NP 

23 


NCA 

80 


NF 

55 


T1 

0.008891 

m 

NX 

6 


NY 

6 


m 

0.01 


CLCA 

0.0052 

kg/m^ 

CLAN 

0.0034 

kg/m^ 

CU 

0.15 


SA 

50000 

m2/kg 

SHO 

0.000044 

-m2 


Length of cell plate in x-direction 

Length of cell plate in y-direction 

Designed current density 

Utilization of O2 in stack 

Utilization of H2 in stack 

Pressure of cooling air 

Operating pressure in stack 

Inlet temperature of process air 

Depth of fuel channel 

Width of fuel channel 

Number of cooling channels 

Thickness of cell (electrode and matrix) 

Inlet temperature of fuel 

Thickness of cell plate 

Number of plates between two cell 

plates 

Depth of process air channel 

Width of process air channel 

Number of cell plates 

Number of process air channels 

Number of fuel channels 

Thickness of cooling plate 

Finite difference number in x-direction 

Finite difference number in y-direction 

Criteria for convergence 

Catalyst loading on cathode side 

Catalyst loading on anode side 

Utilization of catalyst 

Surface area of catalyst 

Cell resistance at 450 K 


Table 5-1 Input data for simulation of dynamic PAFC stack 
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Definition 


Vairiable Dimension Initial Unit 
name value 


ALFA 

0.5 


DKG 

240000 

A/atm 

R 

8.314 

j/g-mole-K 

Z 

2 

g-equivalent 

FCONST 

96500 

C/ g-equivalent 

NC 

36 


KX 

2.5961 

j/m-s-K 

KY 

51.92205 

j/m-s-K 

TKC 

403.3 

K 

WCW 

0.00559 

m 

WCD 

0.00559 

m 

Y1H2 

0.76 


Y1C02 

0.24 


Y202 

0.208 


Y2N2 

0.782 


Y2H20 

0.01 


RHOP 

2611.01 

kg/m3 

RHOC 

2162.49 

kg/m^ 


1.0467 

kj/kg-K 

ccc 

0.841547 

kj/kg-K 


Transfer coefficient 

Constant to calculate limiting current 

density 

Cas constant 

Ntimber of Faraday equivalents transferred 
Faraday constant 

Ratio of cooling air to air cons\mied in 
stack 

Effective thermal conductivity in 
stacking direction 

Effective thermal conductivity on the 
cell plate 

Inlet cooling air temperature 
Width of cooling channel 
Depth of cooling channel 
Mole fraction of H2 in anode inlet 
Mole fraction of CO2 in anode inlet 
Mole fraction of O2 in cathode inlet 
Mole fraction of N2 in cathode inlet 
Mole fraction of H2O in cathode inlet 
Density of cell plate 
Density of cooling plate 
Eeat capacity of cell plate 
Heat capacity of cooling plate 


Table 5-1 continued 
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It is not necessarily true to say that the final steady state is 


reached at 280 seconds. Actually at that time the calculated 
temperature of the stack is still fluctuating. This fluctuation is 
caused in part by the damping characteristics of the system and in part 
by the computation errors. From Figure 5-15, the final steady state may 
be considered to be reached at 43 seconds when the mean temperature 
reaches a steady value. 


s"'- 
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CHAPTER 6 


TRANSIENT RESPONSE OF THE PAFC SYSTEM IN LOAD CHW^GING PERIODS 

The objective of this chapter is to combine all of the programs 
developed in the previous chapters to determine the performance of the 
power plant and its components in the following transient conditions: 
sudden and gradual change from full power load to partial load, from 
standby to partial load, and from the shutdown to standby condition. 

The final computer program was then used to simulate the 
performance of a 7.5 megawatt power plant (Figure 1-4) subjected to a 
power reduction linearly ramping from full load to partial load. The 
results obtained include: the temperature, pressure, and flow rate of 
the components at various locations in the system; the transient 
temperature and current density distributions in the fuel-cell stack; 
the dead time for major components considered in the simulation; and the 
amount of heat and electric energy output. 
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6.1 i^ssumptions 

For simplicity, several assumptions were made to simulate the load 
change responses. The five assumptions made for the fuel cell stack 
subsystem in Chapter 5, are also applicable to the simulation of the 
whole system. In addition, the following assumptions were also made for 
the transient state study. 

1. There are voltage and power constraints, resulting from the fact that 
the platinum and carbon components can corrode at cell potentials above 
0.8 volts. This maximum allowable voltage determines the minimum load 
at which the system can operate. 25% of the fuel load has been chosen 
as the minimum load. 

2. The loads linearly decrease or increase with time. 

3. The conversion rate of methane in the reformer is fixed, as well as 
the carbon dioxide to carbon monoxide ratio. Therefore, the amount of 
input methane is proporitional to the amount of hydrogen needed in the 
stack. The temperature level will be monitored to determine the flow 
rate of methane into the burner. 

<1. The inlet air which flows through the air blower, is proportional to 
the inlet fuel at the burner with the excess ratio fixed. 

5. The steam to carbon ratio is fixed and transient changes in the 
economizer and steam generator are not considered, therefore the inlet 
steam is proportional to inlet methane. The transport lag of water 
pumping and steam generating is neglected. 

6* The recuperator, ZnO bed, and cooling power are not considered in 
this transient study. 

Transport lag is not considered in the piping which combines 
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6.2 Program Description 


A computer program for simulating the transient behavior of the 
PAFC power plant was developed. The steady state conditions for the 
system presented in Chapter 2 were used for the initial conditions in 
the transient state calculation. Efficient, i.e. numerically stable and 
adequately accurate computer models with relatively short computation 
times, were utilized for the individual components and were described in 
Chapters 4 and 5. Finally, by means of the "on-off" flags, which 
controls computations between components, simulation of real-time 


computation can be done on a batch-type computer. 


This Fortran program, which includes 415 subroutines and 20 
functions (descriptions of these subroutines and functions are in 
Appendix 8), was used to simulate the system's transient responses, with 
the required power load being decreased from full load to partial load. 


Although no experimental data was available, this simulation 
provides information concerning the optimal design, operation 
requirements, and control procedure within the load changing period. 


The main (executive) program calls three subroutines which are 
DATAIN, STHAIN, and STTRAN for readii^ input data, calculating initial 
steady state performances, and simulating transient state responses, 
respectively. The flow chart of main program is shown in Figure 6-1. 


The steady state solution was discussed in Chapter 2. The 







START 




READ input 



Figure 6-1 Flow chart of Main program 
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following section will concentrate on the transient responses. 

To decrease the computation time, several assumptions have been 
made. They are 

1. neglet the temperature difference between the stream and the catalyst 
in the reactors. This was discussed in Section <1.41; 

2. decrease the finite-difference section. There are two advantages 
resulting from this modification: first, fewer equations must be solved; 
second, because the maximum stable time interval is decided by the time 
needed to travel one section (except in the reactors). A decrease in 
section number leads to an increase in the maximum stable time, and then 
an increase in the time interval. However, this procedure will also 
lose some accuracy in the final results; 

3. increase the time interval in each component to its maximum stable 
value: If there are two streams in a component (like the tube-side and 
shell-side streams), the greatest common divisor (GCD) of the two 
maximum values is chosen to be the fixed time interval. Thus, the 
system time interval is estimated from the GCD of these components' time 
interval; 

4. neglect the transport lag in the pipes connecting components and in 
the auxiliary components. 


“PAGE MISSING FROM AVAILABLE VERSION’’ 
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g.3.1 Simulation of the Real-time Procedure PRECEWW6 J^A^E BLANK NOT FK.RSED 

To simulate the real-time computation on a "batch" computer, the 
program uses controlled clock time and "on-off" flags in the computer 
coding. There are four flags which will be discussed as follows: 

(1) . IFLRG: flag to indicate the load changing period, 0 indicates 
outlet power under ramping, whereas a 1 indicates that the required load 
has been reached. 

(2) . IFLG: flag for the registering of one of the components so that the 
component may be evaluated for its transient response. Uhere a 0 means 
no registered reaction of the component and a 1 indicates registration. 

There are two situations for ”0" being assigned to a component, they are 
the initial steady state and the final steady state. 

(3) . NFLG: NFLG is used when the calculation will be "operated" at some 
system processing time. Because the time intervals for calculation of 
the components are multiples of system^ s time intervals, at some 
specific system time, some components may not be calculated even if the 
IFLG of these components are on. 

(4) . MFLG: flag is used to indicate the condition of the outlet (s) in a 
component. There are three possible values, 2 for initial steady state, 

1 for transient state, and 0 for the final steady state. By means of 
these "on-off" flags, the calculations among these components 
(subroutines) can be controlled, and the results printed at specific 
transient times to present the components and system conditions. 

Uhen the system time-0, all tlFLGs are equal to 2, IFLGs are equal 
to 0 except for the IFLGs of Heat Exchanger-2, Air Heater, and Fuel cell 
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stack (see Figure 4-1), NFLGs are equal to 0 , and IFLAG is 0 also. The 
MFLGs will be changed to 1 if the corresponding outlet's temperature is 
greater (or less) than the initial steady state value determined by 
CRT2. CRT2 may be some specific value of a temperature criteria or some 
percentage criteria. The MFLGs will be changed to 0 if a change of the 
corresponding outlet temperature is smaller than Y»CRT1, where Y is the 
outlet temperature and CRTl is another convergence criteria. This 
situation will last for some period of time (5 seconds in the example). 
Then a new IFLGs' value will be assigned according to the previous 
component's MFLGs. In addition, the NFLGs' value will depend upon the 
corresponding IFLG and component's time interval. 


6.3.2 Transport Lag Consideration 


Because this is a gas system, the transport lag of streams in 
pipes, reactors, and channels cause a "time delay". A typical example 
would be one pipe (finite-difference) section (in the heat exchanger for 


: t 

I 


instance)^ transporting gas from one location to another. An input flow 
rate change will not immediately transport to the exit. This situation 




can be characterized by specifying the inlet stream as I and the outlet 
stream as J. The flow rate in J will be "delayed" from the flow rate in 
stream I by the time to travel down the section length 1. Such a dead 


mi 


time can be simulated on a computer by a "time delay" procedure. This 
can be constructed in several ways, the most direct of which would be 


the "bucket brigade" approach, a self-explanatory .term (Ref. 31). 


The operating principle of this approach. is to allocate NTN spaces 
to the delay channel, where NTN=time delay/(IDH x SDT) where IDH is a 


1 


multiple of the system time interval (SDT) for a specific component. 
Instead of feeding the input signal (flow rate for instance) in at the 
front end, and moving the values in each space one position toward the 
exit end, and finally reading the value in the last space as the exit 
value ("bucket brigade"), the technique does not move the values but the 
pointer for the readin and readout moves from space to space. 


Each difference section is treated the same way. For simplicity, 
the mean flow velocity is used to estimate the transport time (time 
delay) of one section, then the WIN for every section is the same. 
However, there may be two different values of NTH for one component, one 
is of the tube-side gas, and the other is the shell-side gas. 


6.3.3 Determination of Final Steady State 


Determination of the final steady state is a problem because of the 
fluctuation of the calculated values. However, there are several 
criteria to use, 

(1) . calculate the final steady state performance using the required 
input data, then when transient responses reach 99.5% or more of these 
values, the solutions are considered final steady state solutions. 

(2) . decide by the temperature damping ratio - if the ratio is less than 
0.1%, e.g., the final s.s. is assumed to be reached. 

(3) . decide by the temperature change between the continuous time period 
~ if the change is less then 0.05**K, say, and lasts for a period of time 

5 seconds), it is assumed be the final steady state. 

F’or simplicity in the computer programming, the last criterion was 
selected for this simulation. 




6.3.4 Flowchart of the Progrram 


The flowchart of the executive program is given in Figure 6-2, 
which shows how the "on-off" flags function and how the program executes 
the computations for the transient period. 



Figure 6-2 Flow chart, of PAPC system transient state si m ulation 
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6.^ Example and Results 

The example discussed in this section involes the simulation of the 
PAFC system response for an output power decrease from full load to 90% 
load. The assumed ramping rate is 5% per second for this power 
decreasing operation. Appendices 6 and 7 provide the definitions and 
values of the input data for the steady state and transient state 
analyses, respectively. Some of the input data used here for the steady 
state simulation differs from that of the previous simulation in Chapter 
2. The values are listed in Appendix 6.2. All of the input data is in 
SI or SI acceptable units. 

This sample example was run on the NASA-LeRC IBn/370 computer. It 
required 115.39 minutes of CPU time to reach the final steady state 
which occured at a system time of 14.38 minutes. The criteria used for 
determining the final steady state was discussed in Section 6.3.3. 



iji ’ 

: 



The results include the operating conditions in each 
finite-difference grid (in the fuel cell) or section (in other 
components) at specific transient times, and initial and final steady 
state solutions. Portions of the solutions for the transient time and 
both of the steady state periods are in Appendix 9. Some of the results 
will be discussed here. 

Figure 6.3 shows .the temperatures of some important locations with 
respect to time, the locations are: the input fuel after heat 
exchanger-1 (A8); the mixture of steam and fuel (AlO); the reacting gas 









202 



j 

I 

i. 


outlet of reformer (A12); the outlet of the high temperature shift 
converter (M5); the shell-side outlet of the heat exchanger (M6); the 
outlet of the low temperature shift converter (A17); the shell-side 
outlet of the air heater; and finally, the anode outlet (Bl). The 
nomenclature for the streams refers to Figure 1-<I. Note that most of 
the stream temperatures fluctuate somewhat before the first minute. 
Because there are several process loops in the system, the temperature 
at any location is affected by other changes in the system. This is 
specially true for the initial changes which result from adjustments in 
the flow rates at four Inputs; air for combustion and the fuel cell 
stack, and fuel for the reformer and the fuel cell stack. 



Table 6-1 shows the operating conditions of the fuel cell stack in 
a transient period, the conditions include the operating voltage, mean 
current density, average temperature in the stack, anode side outlet 
temperature, operating pressure, and output power. Because a fixed 
anode inlet composition was assumed (Section 5.1), the fuel cell stack 
subsystem is independent from other components in the calculations. 
Therefore, the time needed for the fuel cell stack to reach the final 
steady is much faster than that of the whole system, actually, the fuel 
processing subsystem. The time period is 3.3 minutes compared to 14.38 
minutes for whole system. In Table 6-1, it is evident that the current 
density and operating voltage are changing dramatically in the first 2 
seconds when the manipulations (of input flow rates) occur, whereas no 
significant change in temperature appears after that time. Table 6-2 
also lists the change of voltage and current density in the first two 
seconds. Note that the linear change in these two sets of data has 



? 

I 2.052 0.685 0.2S7 





Operating voltage and mean ctirrent density of PAPC 
stack changing with time - small time scale 
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resulted from the linear manipulation of the inlet flow rates. In 
conclusion, the current density is initially affected by the 
concentrations (flow rates) of fuel and air during the adjustment 
period, and then by the temperature distribution in the stack. The 
resulting current density distribution also affects the temperature 
distribution at the same time. 


It is difficult to see both of the distribution changes in Table 
6-1. The temperature and current density at each finite-difference 
section of the middle strip in the symmetric plate (see Figure 2-17) are 
given in Table 6-3 and 6-4, respectively. Both of the changes at the 
outlet sections are greater than those of the inlet sections, they are 
2.49 to 0.76 *K for the temperature and 0.0443 to 0.0306 A/cm* for the 
current density. 


Appendix 10 provides a listing of the initial and final steady 


m: I 
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state solutions for the current density and the temperature 
distributions. 


The assumption that the same concentrations of hydrogen occurs 
after the reformer and shift converters is not accurate (see the listing 


of the initial and final steady state performances). However there is 


only a 1.69% ((0. 5928-0. 5829)/0. 5829) difference after the reformer and 
a 2.35% difference ((0. 7000-0. 6839)/0. 6839) after the shift converters. 
The same conversion assumption is essential to calculate the inlet fuel 
flow rate for the first step. However, this does not greatly affect the 
anode inlet composition, because the cooling tower before the reservoir 





Time - - sec 


Finite-difference section along air flow direotion 


c.d.- 

A/cm 2 

1 

2 

3 

4 

5 

6 

0 

.3074 

.3204 

.3357 

.3397 

.3431 

.3311 

0.324 

.3034 

.316 

.3307 

.3343 

.3372 

.3252 

0.756 

.2963 

.3082 

.5223 

.3258 

.3285 

.3165 

1.08 

.2903 

.3018 

.3153 

.3184 

.3210 

.3096 

1.296 

.2866 

.2976 

.3108 

.3136 

.3160 

.3047 

1.512 

.2829 

.2936 

.3063 

.3090 

.3112 

.2998 

1.836 

.2773 

.2875 

.2997 

.3020 

.3058 

.2925 

2.052 

.2749 

.2849 

.2969 

.2990 

.3006 

.2893 

4.104 

.2755 

.2855 

.2973 

.2992 

.3004 

.2884 

60.48 

.2763 

.2856 

.2970 

.2990 

.3000 

.2881 

198 

.2768 

.2865 

.2972 

.2985 

.2990 

.2868 


Table 6-4 Cxirrent densities of the middle portions of the PAFC stack symmetric plate 
(see Figure 2 - 1 6) changing with time 
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can provide some adjustment. 

Although there is some CPU time spent in the iteration of the 
steady state solution, there is still too much CPU time spent in 
obtaining the transient solution. Other than running this program on a 
more powerful "supercomputer" (like CRAY-1), it is necessary to improve 
the computation time by adding some simplifying assumptions which will 
be discussed in the following chapter. 

This program can be modified to simulate the startup and shutdown 
responses of PAFC system, but the procedures for these operations should 
be standardized beforehand. 



CHAPTER 7 


CONCLUSIOtl 

7.1 Summary 

1. The fuel cell is an attractive option for electrical energy 
generation due to its high efficiency and lack of disturbance to the 
environment. Because of its flexibility in size and operating 
characteristics, the fuel cell can also be placed on-site at the point 
of end use. In this location, the fuel cell not only povides 
electricity but thermal energy as well. 

2. Phosphoric acid fuel cells (PAFC) are the most advanced of all of the 
fuel cell technologies. Pilot plants up to 4.8-hU have been 
successfully operated. 

3. The PAFC system has three subsystems, which are the fuel processor, 
the fuel cell stack, and the power conditioner. The fuel processor 
converts a hydrocarbon fuel to a hydrogen-rich gas that is fed to the 
fuel cell stack to produce DC power. The power conditioner transforms 
the DC power to AC power compatible with user's requirements. 

4. The accuracy of digital simulation depends upon the mathematical 
mdels, parameter estimation, and numerical methods used. 

5. There are two modeling methods used in this, simulation, one is the 
lumped model and. the.otlier..is the. distributed model. The former uses 
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the thermodynamic balances to set up algebraic equations, whereas the 
latter use kinetic and thermodynamic expressions to form simultaneous 
ordinary-differential equations. The finite-difference method is used 
to solve the differential equations. 

6. In distributed modeling, the reformer is treated as a combustion-gas 
heated heat -exchanger type catalytic reactor; the shift converters are 
treated as adiabatic catalytic reactors; and counter-current double pipe 
heat exchangers are used to model the heat exchangers in the PhFC 
system. 

7. The steady state simulation of the reformer will be limited by the 
uncertainties of the rate expression and estimation of the heat transfer 
coefficient. 

8. The regenerative type reformer recycles the product gas to save 10 to 
15% of the energy in the combustion gas, which can be used to evaporate 
the water for reforming. In addition, this can be accomplished by means 
of an external heater using product gas to evaporate part of the 
required water. 

9. The performance of the PAFC is controlled primarily by the rate of 
oxygen reduction on the platinum. 

10. For the design of fuel . cell systems, modeling is needed to obtain 
the local , current generation .. as a function of position, flow rates, 
temperatures, feed concentrations, voltages, etc. Non-uniform current 
generation results in non-uniform heat generation in the cell. Hore 
^at siwuiid 'be removed at the locations ' of highest cuiVent density and 
thus highest heat generation . This gives rise to various possibilities 
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for devising better cooling schemes for optimum performance. Among the 
proposed cooling schemes in the cooling plate, the branched channel can 
obtain the best cooling effect but suffers from structural deflection 
and complication, however, the varying-width scheme can provide good 
uniformities of both current density and temperature and has a simple 
structure . 
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11. In the analysis of the fuel cell stack, both the concentration 
gradient and cooling scheme determine the uniformity of performance of 
the cell plate. High utilization of hydrogen or oxygen will increase 
the concentration gradient. 

12. The models for steady state performance in the reformer, fuel cell 
stack, and the whole PAFC system agree very well with the experimental 
results. 

13. The developed steady state model can be used in the design of the 
components and system of the PAFC. The model is also used for further 
simulation of the transient state. 

14. A combined algorithmic-heuristic approach was developed to 
synthesize the heat-exchanger network in PAFC system. There was a 10% 
cost improvement from this synthesis. 



15. Load change occurs in the operation of the power plant in both 
residential and commerial applications. Some transient operation 
methods and assumptions were developed in order to simulate the PAFC 
system response under transient conditions. 

16. The dynamic study of each component was the first step in the 
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simulation of PAFC system's transient response. The distributed model 
was used to estimate the response of each component. 

17. Usually, the stable time Interval is a problem in a digital dynamic 
study. A small value will require a large amount of computation time 
and a larger truncation error will result, and the larger time interval 
values will lose information in the initial time interval and may result 
in a divergent solution. 

18. The maximum value of the time interval in one component depends on 
the time required to travel one finite-difference section in a 
non-reacting component, or the sensitivity of the reaction rate to the 
temperature in the reactor. 

19. Neglecting the kinetic considerations in the PAFC dynamic study 
results in no accumulation of components in each section and also in the 
equilibrium calculations in the transient state computation. Because 
the time constant of reaction rate process is less than that of the 
material transport and much less than that of the heat exchar^e, this 
assumption is acceptable. 

20. The length of the computation times prevents the use of a real-time 
computer in the simulation of PAFC system's trnasient responses. By 
neans of the ”on-off" flags and a process controlled clock the real-time 
computation procedure was simulated on a batch-type computer. 

21. Although no experimental data is available to test the transient 
state simulation of the components and system, the accurate steady state 
simulation, and the smooth transient responses, indicate reasonable 
dynamic results. 
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22. The computer prcgrem developed Is useful tor esteblishins the 
controller settings, to design a procedure for pover ramping, and to 
study the transient responses. 
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7.2 DISCUSSION ftt^D RECOffilANDMIONS 

1. One way to decrease the computation time in the reactors is to assume 
no accumulation of components in the finite-difference section. This 
assumption was used in the simulation of the fuel cell stack in the 
transient study. There are two advantages resulting from this 
assumption: one is the simplification of the differential equations by 
deleting the concentration terms, and the other is the increased time 
interval. If the transport lag in the reactors is neglected, further 
increases of time interval are possible. A comparison was made with 
these simplified conditions in the reformer, which assumes a fixed wall 
temperature profile. The operation definitions, assumptions, and stable 
time intervals are shown in Figure 7-1, where Model I was used in 
Chapter 4, Model II was used in Chapter 6, and Model III was proposed 
here with previously discussed assumption(s). Model III-A, taken from 
Chapter 5, also includes the assumption of no accumulating components, 
in addition. Model III-B also assumes no transport lag in the reformer. 
It is evident thiat the stable time interval is dramatically increased 
500 times by using Model III-B instead of Model I, or even in comparison 
with Model II (40 times). A futher reduction of computation time 
results from fewer variables in Model III. However, the results of Mode 
III-B lose the information for the design of the operation procedure 
(ramping rate) and controllers. If we treat the PAFC stack as a 
reactor. Model III-B can save a large amount of computation time. 

2. Because of the iterative trial -and-errpr procedures in the steady 
state calculation, the accuracy of the -final solutions depends on the 
convergence criteria. Loose criteria may result in an unstable 


I Descriptions Model I 


Additional 

Assumptions 


Model III 


(all assuTiption 
in Section 
2.1.3 are 
available) 



difference 
between 
stream and 
catalyst 



1. no tenqjeratura 1 , no temperature difference 


between stream and 
catalyst 

2. no accumulation of 
components in the finite- 
difference section 

3. A: transport lag of each 
section is equal to one 
time interval 

B; no transport lag, time 
interval is greater than 
transport time through 
the reformer 


Time Interval 0.0072 sec. 


0.09 sec. 


Variables 


l.mole fraction! 1 .mole fraction 


of CH. 


of CH/ 


2. reforming gas 2. reforming gas 

density density 

3. reforming gas 3. reforming gas 

tec^p . ten^P . 


0.18 sec. 


1 .reforming 
gas temp . 


3.6 sec. 


1 .reforming 
gas teiijp . 


4. catalyst 



General Descriptions; height of reformer; 40 ft 

2 o 

tenperature profile of wall; 1300+17x-0.2x , P; x; height, ft 
inlet tenperature of reforming gas; 687 *P 
finite-difference sections nvunber; 21 
S.S, outlet conversion; 82?^ 
kinetic expression: 

where T in *R 

• initial condition; input flow rate decrease to ^ 0 % of 

original, rate; ^ 0 % per 0.18 sec. 

^gure 7-1 Description of different models to conpare the stable time interval 
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calculation (divergence) in the transient state calculation which uses 
the steady state solution as the initial conditions. The convergence 
error may accumulate with the initial manipulation changes and result in 
a divergent solution, h preliminary step can be taken to "stabilize" 
the steady state solution. The dynamic program with the steady state 
solutions is run, without any initial changes, until the difference 
between the calculated steady state solutions and the "stabilized" 
solutions is reduced, which means that the accumulation terms approach 
zero. This require some extra computation time, however. It is 
important to start from the "real" steady state conditions (the initial 
conditions) . 

3. Non-uniform temperature distribution in the stack results from 
non-uniform current density distribution. The non-uniform concentration 
profiles cause these non-uniformities. Use of Zee plate or hexagonal 
plate stated in Section 2.2.2 is one solution. The combinition of new 
cell plate configurations, effective cooling shemes, and suitable 
coolant is a goal for future research. 

4. In order to reach a balance between accuracy and computation time, 
this system's dynamic model neglects some small changes which can be 
considered in future work. The steady and fixed transport lag in each 
component can be reconsidered as a function of position and time, the 
transport lags in pipes and auxiliaries can be estimated and added to 
the simulation, k more important consideration is that the conversions 
in reformer and shift converters should not be fixed. ftll of the 
calculation can be based on this data which will be updated with the 


nost recent values. 
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5. Good models can also be used to explore behavior where only limited 
experimental data are available. For example, atmospheric pressure data 
can be used to model fuel cell performance at elevated pressures. 


6. Digital simulation is a powerful tool for the design of the process 
plant, and it may be expected that the day will come when engineers will 
substitute pilot plant establishment with simulation on digital 
computers. 




7. Mthough PAFC power plant operation is not as critical as that of a 
nuclear power plant, it is necessary to establish a simulator which can 
be used for the design of safety features, h number of "what if" 
questions (sometimes problems) can be answered with this type of 
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simulator. The program developed here is really the first step in the 
development of a PAFC power plant simulator. 


8. Other future research could be on (1) the estimation of wall heat 
transfer coefficients in a catalytic reactor, (2) the kinetic expression 
of the demethanation reaction, (3) a new algorithm with decomposition 
procedure for the heat exchanger network, (d) the reliable prediction of 
the current-voltage relationship. 
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l^ppendix 1 


Minimm Steam to Carbon Ratio 


A possibility exists that during reforming elemental carbon may 


form and deposit on the catalyst bed. Primary control is by proper 


selection of a suitable steam/carbon ratio. 


Kany reactions are possible in this system. A few, including the 


more likely are as follow (Ref. A-l-l): 


CH 4 = C + 2Hg 


(A-l-l) 


C + HgO = CO + H 2 


(A-1-2) 


CH 4 + H 2 O = CO + 3H2 


(A-1-3) 


CO + H 2 O = CO 2 + H 2 


(A-1-40 


2C0 = C + CO 2 


(A-1-5) 


' CO 2 = CO + 1/2 O 2 


(A-1-6) 


H20 = H 2 + 1/2 02 


(A-1-7) 


2 CH 4 = C 2 H 6 + H 2 


(A-l- 8 ) 


Consideration of the values - .. of . these equilibrium constants 


indicates that, at normal reformer operating temperature, reactions 


(A-1-6), (A-1-7), and (A-l- 8 ) can proceed only to negligible extents. 


and hence that O 2 and C 2 Hg cannot be appreciably present at equilibrium. 


If no carbon is to appear., in the . equilibrium mixture represented by 


these three reactions, it is necessary to add sufficient steam 


so that the activity ratio' equal to or greater than 

Kl, the ec^ilibrium constant - of 'reaction -(A-l^l-)^'’- and that the activity 


ratio (a^j^-Xag l/a’g' ^ be "equal"!' to ^ ' 0 ^^ thsih K 2 , the equilibrium 


constant 3 of reaction ■‘<A-1-2)V ' If the first ratio is greater than Kl, 
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then carbon added to such a system can react with hydrogen and form 
methane until the ratio is reduced to Kl. If the second ratio is less 
than K2, then carbon added to such a system can react with steam until 
the ratio is elevated to K2. 


Uhen the ratio of steam to methane in the feed is sufficiently high 
so that carbon cannot be present at equilibrium, the equilibrium 
composition may be calculated by considering only reactions (A-1-3) and 
(A-1-4), which involve all the significant reactants in the absence of 
carbon. In order to determine the minimum steam ratio required for 
freedom from carbon, a constrained minimization problem was set up; 
min. x(3) 
subject to: 


1. < 3x<l)+x<2) P 

< l-x(l) ) < l+x(3)+2x(i) > > Kl 


2. ( 3x(l)+x(2) ) ( x(l)-x(2) ) P 

< x(3)-x(l)-x<2) ) < l+x(3)+2x(l) ) < K2 


3. ( x(l)-x(2) ) ( 3x<lHx(2) >» 

< l-x(l) ) < x(3)-x<l)-x(2> ) ( l+x(3)+2x(l) >» = K3 


x(2) (3xg)+x(2) ) 

. _ < x(l)-x(2) ) (x(3)-x(l)-x(2) ) = K<l 


5. x(2) (H^x<3)+2x<l) ) 

(x(l)-x(2) >* P > K5 


6. x(2)-x(l) < 0 


where x(l): moles CH^ converted by reaction (A-1-3)' 
’ x(2): moles CO converted by reaction (A-1-4) 

x(3): steam to carbon ratio 
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Kl, K2,..., K5: equilibrium constant of reaction (A-l-l), 

...» (R-1-5) 

So at equilibrium (assume initial CH is 1 mole), 
moles CH 4 : 1 - x(l) 

x(3) - x(2) - x(l) 

CO; x(l) - x(2) 

CO 2 : x< 2 ) 

H 2 : 3x(l) + x(2) 

Total moles: x(3)+l+2x(l) 

This minimization problem was solved by COUPUTE computer program 
(Ref. A-1-2). The COMPUTE code solves the constrained optimization 
problems using mixed penalty funtion* together with Hooke and Jeeves 
pattern search method for extremization. 

Each equilibrium constant was calculated utilizing a correlation of 
the form 


a b c 

In K = T’ + T* + T + d. (A-1-9) 

This correlation is found to be a very good fit over the range 800-2000 

•! By mixed penalty functions, we mean if the first 1 constraints are 
inequalities and constraints ( 1 + 1 ) to m are equalities, our problem 
becomes: minimize 

^ = f o (x)-K ZZ ln(gi (x)) + 1 /K ZZ <gi (x))* . The function 
^ (x,K) is then minimized for a sequence of monotonically decreasing K 


> 0 . 












Minimum S/C ratio at different temperatures and pressures 
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Appendix 2 

External and Internal Effectiveness Factors 
in Reformer and Shift Converters 

Consider the fixed bed of catalyst particles shown in Figure A-2-1. 
In general the easiest parameters to measure in such a reactor are inlet 
and outlet conditions and perhaps a temperature profile. 



A-2-1 b«d rvMetor 

The surface reaction rate expression explicitly contains 
concentration and temperature terms which must be evaluated at the 
surface. These conditions may be very different from the measurable 
conditions. 

If we let the observed (sometimes termed "global") reaction rate of 
component A in a fixed bed reactor, , be ‘a' function' of the measurable 
quantities, ' 



230 



*obs= ^global"® Rs(C ,C Tq) ET (A-2-1) 

where a = active catalyst surface area/unit 
volume of reactor bed 
(sq.ft, catalyst/ cu.ft. bed) 

R = heterogeneous (surface) reaction rate, 

(#moles A/(hr) (sq.ft, active catalyst)) 

E = internal "effectiveness" factor 
qr = external "effectiveness" factor 
CAo, CTBo,..., To = bed avg. comp, and temp. 

The internal effectiveness factor, E, accounts for the differences 
between concentration and temperature at the external surface of a 
catalyst particle (CAs, Ts ) and on the inside of. the catalyst pores (CA 
, T). 

The external effectiveness factor, T » accounts for the differences 
between conditions at catalyst surface (Cas, Ts ) and the bulk gas (CAo, 
To). 

Internal Effectiveness Factor 

An examination of Figure A-2-2 shows that the concentration of 
reactant A is at maximum at the pore mouth and then decreases down the 
length of the pore. Thus, the surface reaction rate must vary toward 
the center- of the catalyst particle. He need to evaluate the difference 
between the actual ..reaction rate in the pore and the rate when the pore 
®*o.uth concentration prevailed throughout the pore. - 
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figax^a llBdtl of « estaJjtt b«ad vith a alcropor* 


This difference is evaluated in terms of the internal 

"effectiveness factor" which is the ratio of the rate of the process 

occuring in the pore to the chemical reaction rate which would occur if 

the inner pore surface was totally exposed to the pore mouth 
/ 

concentration. 

Thus, 


Rate of process occurring in the pore 

E = Rate of surface reaction 3 gas phase conds. 


External Effectiveness Factor 


No matter how active a catalyst particle is, it can be effective 
only if the reactants can reach the catalytic surface. The transfer of 


reactant from the bulk fluid to the outer surface of the catalyst 
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particle requires a driving force, the concentration difference. 

The concentration of reactant is smaller at the surface than in the bulk 
fluid. Hence the observed rate, the global rate, is less than the 
intrinsic rate evaluated at the concentration of reactant in the bulk 
fluid. The same reasoning suggests that there will be a temperature 
difference between bulk fluid and catalyst surface. External 
effectiveness factor is to evaluate these differences. Thus 

^ _ heterogeneous surface reaction rat e a surface conditio ns 

I ~ heterogeneous surface reaction rate a bulk gas conditions 
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Effectiveness Factors of Catalysts Used in the Reformer and Shift 
Converters 


There are several different kinds of of catalysts used in the 
reformer and shift converters. For instance, in the reformer, there are 
four kinds of catalyst for experimental or industrial use: Girdler's 
G-56B (Ref. 9 and 10), I.C.I. 57-1 (Ref. 11 and 12), Catalyst 100, and 
Haider Tops0e's R13^A/S (Ref. 11). In this study, only one kind of 
catalyst for each kind of reactor is examined, which are Girdler's G-56B 
for reformer, 93% Fe 2 03 and 7% Cr 2 03 (Ref. 13 and 1<J) for high 
temperature shift converter, and Girdler's G-66B (Ref. 15) for low 
temperature shift converter. 


The definitions of variables and algorithm used in the calculations 
of both internal and external effectiveness factors are in Table A-2-1 
and A-2-2, respectively. The data of catalysts, reactions, and the 
reacting gases are listed in Table A-2-3. Some catalyst data is 
estimated from similar kind of catalyst. 


The final results are listed in Table A-2-41. The G-56B catalyst 
shows small internal effectiveness factor, which means the concentration 
of methane and the temperature of the reforming gas are much lower down 
the pore than at the pore mouth. Therefore, the. calculation of global 
reaction rate of demethanation reaction in the reformer should consider 
this difference if G-56B is; used as the "catalyst. The assumption, which 
neglects the internal effectiveness in calculating the reaction rate, is 
generally correct in the simulation of shift converter. Usually when 
the differences between the conditions at catalyst surface (Cas, Ts) and 


Variable Description 


Unit 


DP size (diameter) of catalyst 

APR average pore radius 

PP catalyst porosity 

II shape of catalyst 

IM model used in cal. effective diffusivity 

RMA aver, radius in macro region (lM=2) 

RMI aver, radius in micro region (lM=2) 

EPA void fraction in macro region (lM=2) 

EPI void fraction in micro region (ll^2) 

NA stoi, n\xmber of reactant (IM=1) 

NB stoi. munber of product (U1=1) 

TOR tortuosity factor (lM=l) 

DAE effective diffusivity (lM=3) 

RID internal dia, of reactor 

BWE bed weight 

BD bed density 

BTH bed length 

EPS bed porsity or for sphere (Ref. 1 & 2) 

TOP operating ten^jerature 

P operating pressvire 

X conversion 

PL(i) average flow rate of gas i 

i=1:CH4, 2:00, 5sC02, 4:H20, 5sH2, 6:N2, 

DLC thermal conduotiviiy of catalyst 

PACTR factor for cal, effective diff, (option) 

ACT activation energy 

N order of reaction 


ft 

A 


ft 

lb 

Ib/ft^ 

ft 

R 

ATM 


Ib-mole/hr 

7:02 

BTU/hr-ft-R 


BTU/lb-mole 


Table A-2-1 Definitions of variables for calculating 
effectiveness factors 


m 


(1) Check radial temperature gradient; (Ref. 5) 

(2) Check axial dispersion effect; (Ref. 5) 

(3) Check differential reactor approach; 

(4) Calculate superficial mass flow rate of bulk gas (CO); 

(5) Calculate Reynold Number (RE); (Ref. 6) 

(6) Calculate Chilton and Colburn JD factor (JD); (Ref. 3) 

JD=1.82*RB** (-0.51) RE <350. 

JD=0.99*RE** (-0.41) RE >350. 

(7) Calculate Schmidt Number (SC); (Ref. 4* 6, and 7) 

(8) Calculate mass transfer coefficient (KG); 

(9) Calculate (or input) observed reaction rate (ROBS); 

(10) Calculate characteristic length (AC); 

(11) Calculate Damkohler Number (DA); 


RE <350. 
RE > 350. 


(12) Calculate (or input) effective thermal conductivity (DLA); 
(Ref. 4 and 6) 


(13) Calculate Prandtl Number (PR); (Ref. 6) 

(14) Calculate factor BATA; 

(15) Calculate Weisz Number (Arrhenius Group) (RO); 

(16) Use Newton's method to estimate external effectiveness 
fcictor; 

(17) Calculate surface cone. (CAS); 

(18) Calculate surfane temp. (TS); 


(19) Check the cone, and temp, differences eure significant 
or not; (Ref. 5) 


Table A-2-2 Algorithm to calculate internal and external 


effectiveness factors 
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(20) Calculate internal effective diffusivity (DAE); (Ref, 4) 

( 21 ) Calculate Damkohler Group II (PHl); 

( 22 ) Calculate heat generation function (BAT); (Ref, 5) 

( 23 ) CalcTilate max, internal temp, difference; 

( 24 ) Check plots in (Ref, 8) to estimate internal effectiveness 
factor or check PHI (1,/ABS(N-R0*3AT)) (Ref. 5) to see 
internal effectiveness factor is significant or not ( 5 % 
interval) 


Table A-2-2 continued 
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DATA 


Catalyst G-56B $ 3 % Pe203 T / CrzO -^ G- 66 B 


Reaction steam reforming 

water shift (high) 

(low) 

DP 


0.015625 

0.001837 

0.00177 

APR 

305. 

32.5 

175. 

PP 


0.407 

0.625 

0.56 

II 


2 

1 

1 

IM 


2 

5 

1 

RMA 

993. 



RM 

33.2 



EPA 

0.528 



EPI 

0.0792 



NA 




-1 

MB 




1 

TOR 



4 . 

DAE 


0.02286 


RID 

0.0677 

0.002133 

1.41 

BWE 

0.02866 

0.001433 

668 . 

BD 


80 . 

85.52 

76.78 

BTH ■ 

0.1025 

0.0469 

5.58 

EPS 

0.467 

0.39 

0.51 

TOP 

1640. 

1177.5 

888 . 

P 


1 . 

1 . 

1.5 

ROBS 


6.52482 


X 

[1: 

0.054 


0.78 

FLI 

I 0.02145 

0 . 

0.254 

( 

ri 

1 0.0005968 

0.0004966 

0.57 

1 


O.OOI7O8 

0.0004192 

2.626 

( 


0.0653 

0.000854 

2.29 

( 

p! 

0.001984 

0.00038 

15.27 

( 

6 ] 

0 . 

0 . 

0 . 

( 

7 ) 

0 . 

0 . 

0 . 

DLC 

1 . 

0.5 

7. 

PACTR 0.07 



ACT 

36000. 

43200. 

22860 . 

N 


1 

1 

1 


Table A- 2-3 Input data of examined catalysts 
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G— 66b 


G-56B 939^ Pe205 7% Ct20^ 


CAS/CAO 0.99 

tas/tao 0.99 

External 

Effectiveness 0.8872 

Factor 

Max, Internal 

Temp. Diff. -27.8 

Internal 

Effectiveness 0,3 

Factor 


0.990 

0.9999 

1.000 

1.000 

1.01 

0.9999 

0.87 

0.042 

1.00 

1.00 


Table A-2-4 Results of examined catalysts 
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the bulk gas (Cao, To) are small, then the external effectiveness factor 


is close to 1. 
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{Appendix 3 


Heat Exchanger Network 


A combined algorithmic-heuristic approach to the systematic 
synthesis of heat exchangers, heaters, and coolers is proposed. This 
three-phases algorithm may generate, with great ease and considerable 
rapidity, a near optimal exchanger network. All the sample problems 
ranging in size from 4 to 10 streams were examined. 


Results compare favorably with the presently available techniques 
and allow to handle realistically large problems. 


Introduction 


The area of energy conservation that has been receiving increased 
attention is the improved process heat recovery. Any heat recovered and 
reused in the process not only reduces the amount of fuel consumed but 
also lowers the amount of heat rejected to the cooling system. 


Essentially, the synthesis task consists of finding a feasible 
sequence of heat exchangers in which pairs of streams are matched, such 
that the network is optimal as the total cost is minimum. 


The general techniques that have been developed recently for 
solving this problem include the heuristic approach which is based on 
the use of rules of thumb (for example, H/H combination in Ref. A-3-12), 
and algorithmic methods which often involve some established 
optimization principles (for example, Branch-and Bound in : Ref . A-3-6; 
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illustrate the method. The results are compared with the present 
synthesis methods. Mso, more detail considerations of different 
structures are examined and some recommendations are also given. 


Although this proposed method does not favor all the problems, less 
computation time and relatively easier algorithm still give us good 
understanding in this kind of problems. 


‘ Further extensions of the present approach along with decomposition 

for the synthesis of heat exchanger network are recommended. 

I 

Problem Statement 

i The synthesis problem to be considered has been defined by tlasso 

and Rudd (Ref. A-3-9). In brief, there are M hot streams 
Shi(i=l,2, . . . ,M) to be cooled and N cold streams Scj(j=l, 2, ...,N) to be 
heated. Associated with each stream are its known input temperature Ti,. 
output temperature Tti, and heat capacity flow rate Hi. There are also 
available auxiliary steam heater and water coolers Suk(k=l,2, . . . ,p) 
called utilities. The problem is to create a minimum cost network of 
exchangers, heaters and/or coolers so that the desired output 
temperature of each process stream is reached. Heat transfer area of an 
exchanger of this type is expressed by 
A = Q*log(d2/dl)/(UF*(d2-dl)) 
where dl=Tho-Tci, d2=Thi-Tco. 

& S general, the investment cost for the ith exchanger, heater, and 

cooler, denoted by Cei, Chi, and Cci, respectively, can be correlated to 
their corresponding heat transfer area Aei, Ahi, and Aci by the 



total cost of investsient and utility of the network to be minimized can 
be expressed as 

J=d«(E7a*Aei'»*b+ E’a*Rhi**b+ IIa*A.ci**b)+Z IT uk^Sukl (A-3-1) 

ei , hi ci k 1 

For convenience, the following simplifying assumptions, which have been 
used in most previous studies of synthesis of heat exchanger networks, 
are included: the use of countercurrent shall and tube exchangers with 

minimum allowable approach temperature (JlAiiT), no phase changes of 
process streams, and equal values of the effective heat transfer 
coefficients for all exchangers. The illustrative examples in this 
paper are all taken from the literature, and their specifications of 
process streams and design data are summarized in Table 1 and 2. 



Algorithm 
phase I: 

(1) Calculate the upperbound of heat duty on the level of energy 
recovery. 

<2) Calculate the minimum heating capacity in the heater. 

<3) Select Tct«=Hax(Tcti ) and Th'»=Hax(Thj) corresponding 
to cold stream i and hot stream (include heater) j, 
respectively. 

(^) Consider a exchanger or heater deciding the Th^ 

of which Tct« and Th* are the cold outlet and hot inlet 
temperature respectively. 

(5) Use an appropriate heuristic (maximal heat recovery or 
minimal exchanger, area for example) to determine the 
• quantity of heat transferred in this exchanger, while for 
heater the minimum heating capacity is considered. 






4SP1 

All other 
problems 

steam: 

2 

Pressure (ib/in. abs) 

962.5 

450.0 

Latent heat (3tu/lb) 

656.6 

767.5 

Temperature ( °p) 

540 

456 

Cooling water: 
Temperature 

100 “p 


Heat capacity 

1.0 Btu/lb“P 


Maximum water output 
temperature 

180“p 



Minimum allowable approach temperatures: 
Heat exchanger 20 “P 

Steam heater 25 "P 

V/ater cooler 20 “P 


Overall heat transfer coefficients: 
Heat exchanger 15O B 

Steam heater 200 B 

Water cooler 15O B 

Equipment down time 580 h 


Network cost 
parameters 

Annual rate of return 
Cooling water cost 
Steam cost 


Table 2 


150 Btu/hr ft2 '“P 

200 Btu/hr f t^ 'P 

150 Btu/hr ft^ ®P 

580 hr/yr for 4SP1 and 6SP1 

260 hr/yr for all other problems 


b = 0.6 


a = 550 b = ' 
5 = 0.1 
5 X 10-5 $/ib 
1 X 10-5 $/lb 


Design Data 
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(6) Calculate the new conditions for the exchanged system 
and delete the fulfilled stream. 

<7) Repeat steps (3) to (6) until no further matches can be made. 

(8) Heat any remaining Tci to Tcti using heaters, and cool 
any Thi to Thti with coolers. 

The above procedure will generate only one process, but it is 
nearly always close to the optimum. 

Phase II 

(1) Calculate the lower bound of units by the following rule 
<Ref. A-3-3) 

"The minimum number of units is nearly always one less than 
the number of streams required in the problem". 

(2) Eliminate the unnecessary heaters when the minimum capacity 
of heater criterion is violated. 

(3) Eliminate the unnecessary exchangers by evolutionary rules 
<Ref. A-3-7 and A-3-8). 

(-y) Go back to phase I and recalculate the new structure and 

compare to the original one. If improved result is obtained 
go to Phase III; otherwise, keep the original structure and 
go to Phase III. 

Phase II can be omitted when the criterion (stated in (2)) is not 
violated. It also can skip steps (1) and (3). 

Phase III 

(1) Count the number, NC, of "critical" exchangers (which means 
that the approach temperature is equal to MAAT). 



2^19 

(2) If KC=0, the optimal solution has been obtained; 

if NC=1, Golden-ratio search method is used to find optimal 
MART and the results; 

if NOl, Pattern search method is used to find optimal 
MAAT and the results. 


Results and Comparison 

The sample problems listed in Table 1 are solved usin2 this 
algorithm by IBH/370. These results have been compared with other 
approaches and summarized in Table 3. 

Discussion and Recommendation 

Phase III will improve a lot in Problem <5SP1. The Golden-ratio 
search method is discussed in Ref. A-3-10, and the results are 
summarized in Table 41. However, in other problems the improvement will 
not be apparent (sea Table 5). Whether the optinxal ttfLAT will be larger 
or smaller than 20 ®F is dependent on the temperature differences in the 
"critical" exchanger and properties of streams. In conclusion, 20 is 
a good approach of flAAT and Phase III can be omitted for saving 
computation time. 


Cost of 

43P1 

43P2 

53? 1 

?asso and 
Hudd 



38762 

Lee et al 

13638a 

72400 

(acyclic) 

33762a 

JfcCalliard 

and 

Westerberg 

1 36aea 


38762a 

Pho and 
Lapidus 

1 3638a 


38762 

Hathore 

and 

Powers 

13550 



Ponton 

and 

Donaldson 

1 3550 

237l6d 

45155 

Kelahan 
and Caddy 

o 

ON 


38762 

C’ishida 
et al. 

1 3550 

20353c 

3871 3c 

Limhoff 

and 

Plover 

1 3550 

2l644ad 

I5557ac 

38777a 

This work 

lC533e 

2l92ldf 

15512c 

38762 


65?1 

73? 1 

73P2 ICSPi 

37331 

34376 

2S32S 

35730 

35659 

32152a 

23518b 44150 


35*107 

40625 


44=7 3 

35048 




35010 



43=84 

35010 

35735a 

27560a 

4;=4=a 

35010 

3C414f 

25325 

443oSh 


as the cost is corrected for the different design data used here 
cs not feasible c; splitting 

ds cyolic es I-'-JlI: 1.J9 ? 


fs MAAT 24.47 ? 
hs I3JC. heat recovery- 


r.in. 


40 ? 

exchanger area 


Table A-3-3 Comparison witn Previous Studies 


''lAAT (?) Cost of Problem 43?1 


20 

13590 

15.5 

12746 

9«6 

11655 

5.9 

11059 

3.7 

10781 

1.4 ' 

10536 

Table A- 5-4 Search 

for Optimal /AAT of Problem 4S?i 


A!AAT (?) 

Cost of Problem 73? 1 


20 

30350 


22 - 

30754 


25 

30654 


27 

3C557 


30 

30534 


35 

30453 


40 

30414 



Table A- 3-p ' Search’ for Optimal MAAT 


of Problem ~3P1 
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"Cyclic" structure is better than "acyclic" structure espically in 
the problems where the number of hot streams is quite different from the 
number of cold streams. Some algorithmic methods (Ref. A-3-6 and 
A-3-11) can only obtain acyclic solution and should not be considered in 
solving this kind of problems. The most apparent example is <JSP2» the 
result is a three-fold decrease in annual cost. 

"Splitting" structure was discussed by Linnhoff and Flower (Ref. 
A-3-8) and will not be repeated here. Two different answers of <1SP2 
using splitting structure are listed in Table 3. Our solution uses 
usual optimal technique to obtain optimal splitting ratio and is 
referred to that of Linnhoff and Flower's. Although this procedure is 
not included in our algorithm in which splitting structure can not be 
created, it is worthy to note the difference for this special problem. 
Besides, in order to ensure that its temperature difference at the hot 
end does not, in turn, violate HAAT, the heat capacity flow rate chosen 
for the bypass must not be chosen below a certain threshold value. 

Maximal heat recovery or minimal exchanger area approach for 
calculating the duty of heat exchanger in which Qci is larger than Qhj 
is considered in this work (in Phase I). Though it does not make much 
difference in Problem lOSPl (see Table 3 the last two results), it is a 
good consideration to obtain a better solution. Some constraints and 
criteria are mentioned for choosing these different approaches in our 
, algorithm. , 


This proposed algorithm does not guarantee the "best" solutions 
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(neither do the other approaches), but the less computation tims, 
relatively easier algorithm, and almost optimal solution still make it 
attractive to go further. 

Decomposition method is recommended to combine with this proposed 
algorithm and hopefully will give a better result. An element duty Q is 
chosen by the greatest common divisor of overall duty (Qhi or Qci), then 
we substitute each cold or hot stream by a number of equivalent 
pseudo- streams with each one labelled by a duty Q. These pseudo-streams 
are then solved by the procedures described in this paper. A large 
number of pseudo-streams, which exists in the revised problem, will 
result in longer computation time. 


; f 
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{Appendix <1 

Use Laplace Transform Method to Solve 
Transient Convective Heat Transf er in the Channels 

In the mathematical function of the transient convective heat 
transfer along the channels in the fuel cell stack, only one grid is 
considered. Is 


Tw 


Tj-1 -H 



pAoC “I = hs <Tw(0) - T) - lii 


mC 




vhere T = T <x,0) 

B.C. X a 0 T = Tj (0) 

I.C. 0 a 0 T = T' 



where ^ ; molar density 

Ao: cross area per channel 
e : time 

* ; . previous value 

then £T ^ _hs, , ©T 

»0 AqC 5x 


AoC 


Let 


hs _ m 
AoC^ ■ ■ . ^ 


mC 

AoC 


= B, and E= x/B 


(A-4-2) 


then ~ = A <Tw(0) - T> B 21 . 

■80 - A ; Tw (8-xs> d : 


Taking Laplace Transforms, ^ have. If T(x,s) =oCtT<x,t)J, 

-:'.f e <J) 

ST - T(x,o> =A • Tw(s) - aY -Bp 

•'1 " ; dx 

dT ... m A-Tw(S) T1 


<A-fl-3) 


(A-q-4> 


(A-4-5) 
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T = (ATw(s) + T') -1- + Clexp (- ^ x> 

A+S B 

at X = 0 T = Tj-1 (S) 

then Cl = Tj-1 (s) - (A Tw(s) + T') ( ) 

A+S 


(A-^-6) 


T 


TwCa)«A T’ „ / V ( A+S v Tw(s)*A / A+S x 

= -sir'— — lis^ (- — =') 

(I) (II) (III) (IV) 


T’ / A+S X 

A^ exp ^ x) 


(V) 


The inverse Laplace transform of 

(I) = / ® A exp <-Aa) Tw (0-u) du by convolution property 

(II) = T' exp (-A0) 

(III) = exp (- A Ex) Tj- (8- E ) U (8- E ) by 

second shifting property 

where U (0- E ) is the Heaviside unit step function 
U (0- E ) = f 1 when 8 > E 
0 when 8 < E 

(IV) = -A exp (- AEx) f® exp (-A(u- E)) U (u- E) Tw (8-u) du 


(A-4-7) 


'9 

-A / exp (-Au) Tw (8-u) du 0 > E 

0 

0 8 < E 

(@= -T' exp (- AEx) exp (-A(0- E)) U (8- E) 

= -T' exp (-A0) U (8- E) 

. at X = A X , , . 

let D = A x/B 


fTj' exp ( 


-A0) + A I « 

' yo 


exp ,(-Au) Tw (0-u) du 


if 0 <.D 




ro 
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exp (- Mx) Tj (e- D ) +A 

■'0 

If 0 > D 

V. 

if we use trapezoidal formula <2 
approximate the integration, then 


exp (-Au) Tw (9-u) du 


points for 0 being 


small) 


to 


Tj' exp (-A0) + l/2(ft0) (Tw(0) + Tw(0) exp <-A8) 

f or 0 < D 

Tj - 

Tj-l (0- D > exp (- D ) + l/2( D > (Tw(0) + exp ( -D ) 
^ Tw (0- D >) for 0 > D 

if we let H ( = A0) ■!== D but larger than D, 
then at 0=H Tj(H) Tj-| (0) exp (- D > + l/2( D ) 


<A-^-9) 


(Tw<H) + Tw (0) exp(- D )) 

0=2H Tj(2H>^ Tj-i (H) exp ( -D > + l/2( D ) 
<Tw(2H) + Tw(H) exp < -D ) 

0=3H Tj(3H) ^ Tj-i <2H) exp < -D ) + l/2< -D ) 
(Tw<3H) + Tw(2H) exp ( -D >> 


etc. 

for j=l Tj-t is constant and equal to inlet temperature. 

The problem is that B is quite different between process air and cooling 
air, therefore H can not be estimated. That is 

' D1 1 

B = and < D ) cooling « 43^ ( D ) process, 

so if we estimate H from D process then H » D coolings i: 

Therefore, for cooling air. 


(A-^-11) 
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at 0=H Tj(H) Tj (H) exp ( -D ) + l/2( D ) Tw(H> (l+exp( -D )) 


0=2H Tj(2H) Tj (2H) cxp( -D )+l/2( D ) Tw(2H) (l+exp( -D )) 


etc. 

But, if D is close to or larger than 1, which happens in the process 
air channel when we consider larger grid (ax is larger), then 2 points 
trapezoidal formula approximation does not work. 


Reconsider the mathematical model and take Tw as a constant (can be 

the average of the time period), then Laplace transform of Tw term is Tw 

S 

and 


T 


TwA T' 

(t«)(A+S) " A+S 

(I) (II) 


Tj-l(s) exp (- x) 

(III) 


Tw.A 

(S)(A+S) 


/ A+S X T' , 
exp (- — x) - exp (- 


(IV) 


(V) 


A+S 

B 


x) 


(A-4-12) 


The inverse Laplace transforms are the same except (I) and (IV) are (Tw 
<l-exp(-R0)) and <Tw <exp<-R8) -exp ( -RE ))) U (8- E), respectively. 
Then the final solutions are 

r 


Tj 


Tj' exp (-A0)+Twj (1-exp (-A0)) if 0 < D 

, Tj-i (0- D )exp(- AD )+Twj(l-exp(- AD )) 


(A-4-13) 


if 0 > D 


After going through the same treatments as before, the solutions are as 


follows 
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for process air; 

e=H Tj(H) = Tj-i(0)exp( -AD )+Twj(l-exp( -AD )) 
0=2H Tj(2H) = Tj-i (H)exp< -AD )+Twj(l-exp( -AD )) 


where Twj = (Tw(0)+Tw(e+H) )/2. 


for cooling air: 

e=M Tj(H) = Tj (H)exp( -AD )+TwJ(l-exp( -AD )) 
0=2H Tj(2H) = Tj (2H)exp( -AD )+Twja-exp( -AD >) 


Appendix 5 

Input Bata for Simlation of CSO PAFC System Steady 


Variable 

name 

Bimens ion 

State 

Initial 

value 

Performance 

Unit 

Befinition 

TOPPC 


443 

•k 

Operating ten^erature in fuel cell 

UT 


0.8 


Utilization of H 2 in stack 

CB 


600 

mA/cm^ 

Besigned current density 

BNSM 

1 

1216. 

g-mole/hr 

Input mole flow rate of CH4 


2 

0 

g^mole/hr 

Input mole flow rate of O2 


3 

1.36 

g-mole/hr 

Input mole flow rate of CO 


4 

21.8 

g-mole/hr 

Bmput mole flow rate of CO2 


5 

166. 

g-mole/hr 

Input mole flow rate of H2 


6 

0 

l^mole/hr 

Input mole flow rate of H2O 


7 

0 

g-mole/hr 

Input mole flow rate of N2 


TAT 

298 

*2 

PAT 

1 

atm 

3 MRA 

3.0 


POPE 

5.0 

atm 

IPUEL 

1 


EE 2 R 

0.01 


IP 

2 


I 

7 


IKT 

100 


SXA 

100 


ZE 

2.438 

ffl 


Ambient temperature 
Ambient pressure 
Steam to carbon ratio 
Operating pressure of reformer 
Input fuel is methane 
«2 Input fuel is methanol 
•3 Input fuel is naphtha 
Criterion of convergence in system 
trial and error procedure 
•1 adiabatic operation in shift 
converters 

•2 Isothermal operation in shift 
converters 

Number of con^onents in whole system 
Sztra percentage of needed air in 
burner 

Sxtra per eentage of air in fuel 

cell stack 

Bei^t of reformer 
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Variable 

name 

Dimension 

Initial 

value 

Unit 

DX1 


0 

m 

DX2 


0.0457 

m 

DX3 


0.0509 

m 

KO 


10400 


£A 


83736 

J/mol 

HHOB 


1281.477 

Kg/m3 

EPS 


0.487 


S 


0.0762 

m 

DP 


0.001 

m 

DZZ 


0.0762 

m 

CN 


1.3 

m2-*E 

U 


56.78263 

j/m2-s-* 

HA 

7 

0.2 

m2 


10 

0.2 

m2 

KPH 


2 



NRE 

BASPC 

ODTH 

PITCH 

CLH 

H)SH 

IDTH 

PLOAH 

SUEPC 

CLENE 


5 

0.3048 

0.019 

0.0254 

0.00634 

0.254 

0.0142 

0.00016 

0.0136 

0.61 


m 

m 

m 

m 

m 

m 

bj2 

m 


Definition 

Outside diameter of regenerative 
tube 

Inside diameter of reforming tube 
Outside diameter of reforming 
tube 

Hate constant of demethanation 
reaction 

Activity energy of demethanation 
reaction 

Density of packing in reformer 

Void fraction in reformer 

Width of clmbustion gas si^uare duct 

Diameter of catalyst in reformer 

Eei^it of finite— difference section 

Q X A/C min in heat exchanger 

Overall heat transfer coefficient 

in heat exchanger 

Transfer area in E-7 

Transfer area in £>-10 

Number of tube passes 

Number of tube rows 

Buffle space 

0. D. of tube 

Pitch of heat exchanger 
Clearance in heat exchanger 
IJ). of shell 

1. D. of tube 

Plow area in heat exchanger 
Surface area per line 
Length of tube 
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Definition 


variable Dimension Initial Unit 
name value 


SITSZ 


0.5 


Ratio of total inside tube cross- 
sectional area per pass to header 
cross-sectional area per pass 

DTH 


0.7 


Fraction of AT over inlet gas film 
in heat exchanger 

DPD 

1 

0,36 

m 

Diameter of shift converters 

AHRU 

1 

0.66 


Void fraction in shift converters 

APPD 

1 

6.41 


Total surface area of packing 

CLEPD 

1 

1.8 

s 

Length of shift converters 

NTPD 

1 



Number of tubes in shift converters 

NT^ 


140 


Number of fuel flow channels in stack 

PUXiA 


0.433 

m 

Length of fuel channel 

WIDA? 


0.00297 

m 

Width of square fuel channel 

NPHI 


3365 


Number of cell plates 

NTAA 


40 


Number of process air flow channels 

A^tXi 


0.3048 

m 

Length of air channel 

WIDAA 


0^00157 

m 

Width of square process air channel 

SRO 


0.000044 

.n.-a2 

Cell resistance at 450 

SA 


40000 

m2/Kg 

Surface area of catalyst 

CD 


0.15 


Utilization of catalyst 

CL 


0.0075 

Kg/a2 

Catalyst loading 

ALFA 


0.50 


Transfer coefficient 

SN 


2 

^equivalent 

Number of fUraday equivalents 





transferred 

PCONST 


96500 C/g-equivalent 

Faraday constant 

DKC 


240000 

A/atm 

Constant to calculate limiting 


current density 
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Appendix 6 

A6.1 Input Data for Sinnilatlon of Westlnghouse PAPC System S.S. Performance 

Variable Dimension Initial Unit Definition 

name value 


XN 

0.35433 

m 

IN 

0,3048 

m 

XDNSO 

3250 

A/m^ 

UTA 

0.5 


UTH 

0.8 


POPC 

3.4 

atm 

POP 

3.4 

atm 

TKA 

450 

Z 

WFD 

0.001016 

m 

wpv 

0.003046 

m 

NCC 

36 


WE 

0.001016 

m 

T 

0.003302 

m 

NK 

5 


WAD 

0.001016 

m 

WAW 

0.003048 

m 

HP 

32400 


NCA 

80 


NP 

55 


T1 

0.008891 

m 

NX 

6 


NY 

6 


EB 

0.005 


CLCA 

0.0075 

kg/m^ 

CLAN 

0.0034 

kg/m^ 

CU 

0.35 


SA 

50000 


S£0 

0.000044 

_fi--m^ 

ALFA 

0.55 


DEC 

^0000 

A/atm 


Length of cell plate in x-direction 
Length of cell plate in y-direction 
Designed c u rre n t density 
Utilization of 02 in stack 
Utilization of E2 in stack 
Pressure of cooling air 
Operating pressure in stack 
Inlet ten^erature of process air 
Depth of fuel channel 
Width of fuel channel 
Ihxnber of cooling channels 
Ihiekness of eell(elsctrode and matrix) 
Thidmess of cell plate 
Number of plates between two cooling plates 
Depth of process air channel 
Width of process air channel 
Number of cell plates 
Number of process air chazinels 
Number of fuel channels 
Ibickness of cooling plate 
Finite difference number in x-Klirection 
Finite difference number in y^irection 
Criteria for convergence 
Catalyst loading on cathode side 
Catalyst loading on anode side 
Utilization of catalyst 
Surface area of catalyst 
Cell resistance at 450 Z 
transfer coefficient 
Constant to calculate limiting current 
density 






Variable Bisens ion 
nase 

. Initial 
value 

TInit 

Beflniticn 

B 


8.314 

j/(ff-mol)(K) 

Gas constant 

Z 


2 

^-equivalent 

Ihnnber of Baradaor equivalents transferred 

BCCBIST 

NC 


96500 

29 

C/^-equivalent 

Baraday constant 

Batio of cooling air to air consuaed in 
stack 

IX 


1.730735 

J/m.a.I 

Effective thermal eonduetivity in stacking 
direction 

Ef 


51.92205 

j/s.s.K 

Effective thermal eonductivi'fy on the 
cell plate 

TKC 

NTREED 


403.3 

3 

I 

Inlet cooling air teoqierature 
■0 Strai^t or varying width cooling 
configuration 

Branched cooling configuration 

vnuRE 

1 

0.1524 

m 

Piret section length of branched cooling 
e onf igoratian 


2 

0.1016 

m 

Second section length of branched cooling 
configuration 


3 

O. 05 O 8 

m 

Third section length of branched cooling 
configuration 

WCW 

1 

0.00508 

a 

First width of branched or varying width 
cooling configuration 


2 

0.00254 

a 

Second width of branched or varying width 
cooling configuration 


3 

0.00127 

a 

Third width of branched or varying width 
cooling eonfiguratian 

¥CL 

NVABY 


0,00508 

0 

a 

Beptb of cooling channel 
■0 Strai^t or branched cooling 
configuration 

•1 Varying width cooling configuration 

IVA 

1 

2 

0 

0 


First width section nxintber 
Second width section cumber 

Y202 

3 

0 

0.208 


Third width section number 

Hoi* fraction of 02 in cathode inlet 


4 



Definition 


Variable Dinenaion Diitial Ibiit 
name value 


Y2N2 

Y2H20 

BHQP 


0.782 

0.01 

2611.01 

kg/m^ 

HBOC 


2162.49 

kg/m^ 

CCP 


1.0467 

kJ/(kg.Z) 

CCC 


0.341547 

kJ/(kg.K) 

fiAWRA 

1 

964.5670 



2 

2411.417 

m^/m^ 

CPC 

1 

0.879228 

kJ/(kg.Z) 


2 

0.753624 

kJ/(kg.K) 

HHOB 

1 

1281,501 

kg catalyst/m' 


2 

1229.920 

bed 

kg catalyst/ffl 

DP 

1 

0.003048 

bed 

m 


2 

0,001219 

m 

HFS 

2HS 

1 

0.469 

1.8266 

m 


2 

O.6O96 

m 

D1S 

1 

0.09144 

B 


2 

0.0762 

m 

NTS 

1 

92 



2 

92 

- 


Mole fraction of in cathode inlet 

Hole fraction of H 2 O in cathode inlet 

Density of cell plate 

Density of cooling plate 

Heat capaci-ty of cell plate 

Heat capacity of cooling plate 

Specific surface area of packing in high 

tes^erature shift converter 

Specific surface area of packing in low 

teng>erature shift converter 

Heat capacity of catalyst in high 

te a pe ra ttge shift converter 

Heat capaci'^ of catalyst in low 

tes^erature shift cinverter 

Density of packing in hi^ ten^erature 

shift converter 

Density of packing in low temperature 
shift converter 

Diameter of eatsilyst in hi^ temperature 
shift c cover ter 

Diameter of catalyst in low tenperature 
shift converter 

Void Araction in shift converters 
Height of high temperature shift converter 
Height of low tesperature shift converter 
Diameter of high tesperature shift 
converter 

Diameter of low tesperature shift 
co n verter 

Number of tubes in h-t gh tesperature 
shift converter 

Number of tubes in low tesperature 
shift converter 
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Definition 


Variable Dimension Initial IMit 
name value 


16 1 

2 

ZE 1 

2 

3 

4 

5 

D1 1 

2 

3 

4 

5 

D2 1 

2 

3 

4 

5 

D3 1 

2 

3 

4 

5 


6 


6 

1.8268 

m 

4.2672 

m 

2.7432 

m 

5.7912 

m 

5.7912 

m 

0.0452628 

m 

0.0432626 

m 

0.0452628 

m 

0.0190729 

m 

0.0190729 

m 

0.0500101 

m 

O.O5O8IOI 

m 

0,0500101 

m 

0.0253090 

m 

0.0253098 

m 

0.073152 

m 

0.073152 

m 

0.073152 

ffl 

0.0370332 

m 

0.0370332 

m 


Number of finite difference sections in 

hl^ tenperature shift converter 

Number of finite difference sections in 

low tes^eratore shift converter 

Length of air heater 

Length of anode exhaust heat exchanger 

Length of preparator 

Length of heat ezchanger-1 

Length of heat exchanger-2 

Inside diameter of tube in air beater 

Inside diameter of tube in anode 

exhaust heat exchanger 

Diside diameter of tube in preparator 

Inside diameter of tube in heat 

ere hanger-1 

Inside diameter of tube in heat 
exehanger-2 

Outside diameter of tube in air beater 
Outside diameter of tube in anode 
exhaust beat exchanger 
outside diameter of tube in preparator 
Outside diameter of tube in heat 
exchanger-1 

Outside diameter of tube in heat 
exchanger-2 

Inside diameter of shell in air heater 
Inside diameter' of shell in axx>de 
exhaust heat exchanger 
Inside diameter of shell in preparator 
Inside diameter of shell in heat 
exehanger-1 

Inside diasieter of shell in heat 
exchanger-2 


266 


Variable Dlstensicsi Initial 
nane Talue 


SefinitioD 


Nuaber of tubes in air beater 
Ruaber of tubes in anode ezbaust heat 
ere hanger 

Nuaber of tubes in preparator 
Buaber of tubes in heat exchanger- 1 
Nuktber of tubes in heat exchanger-2 
Nuaber of finite difference sections in 
air heater 

Nuaber of finite-difference sections in 

anode exhaust heater 

Nuaber of finite -difference sections in 

preparator 

Nuaber of finite-difference sections in 
heat exchanger- 1 

Nuaber of finite -difference sections in 
heat exchanger-2 


RflW 

8027. 17 

kg/a^ 

Density of vail in heat exchangers 

CPW 

0.502416 

kJ/(kg.E) 

Heat capacity in heat exchangers 

THW 

20.76882 

J/m.s.E 

Theraal ccnd'octiTi'^ in heat exchangers 

ZHR 

12.192 

a 

Height of reformer 

DIB 

0 

a 

Outside diameter of regenerative tube 

D2B 

0.12701 

a 

Inside diameter of reforming tube 

LJB 

0.149992 

a 

Outside diameter of reforming tube 

EO 

10400 


Bate constant of deaethanation reaction 

EA 

83736 

J/aol 

Actlvi'^ energy of demethanatlon reaction 

BHCB 

1281.477 

kg/a^ 

Density of packing in reformer 

EPSB 

0.45 


Void fraction in reformer 

S 

0.2286 

a 

Width of combustion gas square duct 

DPB 

0.004572 

a 

Diameter of cat«Qyst in reformer 

MB 

21 


Number of finite-difference sections in 



kg/a^ 

reformer 

BEWB 

8027. 1725 

Density of vail in reformer 
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Variable Dimension 

Initial 

value 

Quit 

CPWB 


0.6112728 

kJ/(kg.K) 

SABEAB 


669.29133 


CPCR 


1 .025766 

lcJ/(kg.K) 

irate 


0 


HTR 


35 


CTAH 

1 

0.95 



2 

0 



3 

0 



4 

0 



5 

0 



6 

0.05 



7 

0 


TAT 


300 

K 

PAT 


1 

atm 

SMRA 


2.5 


linJEL 


1 


IP 


1 


I 


7 


EXT 


10 


EBR 


0.0025 


U 


56.78359 

j/(m^.s.K) 

HA 


371.04 ; 

«2 



0.971 


PPI 


3.4 

atm 


Definition 

Beat eapacitjr of vail in refozmer 
Specific surface area of catalyst in 
reforoer 

Beat capacliy of catalyst in refonser 
■0 Pirst order reaction kinetic of 
desietbanation reaction 
a>1 Vestin^iouBe kinetic model of 
danethanatiaa reaction 
Humber of tubes in reformer 
Mole fkaetlon of in input fuel 
Hole fraction of CO in input fuel 
Mole fraction of CC^ in input fuel 
Mole iksctlon of B20 in input fuel 
Mole Aractlon of B^ ^ input fuel 
Mole fraction of N2 in input fuel 
Hole fraction of O2 in input fuel 
Ambient tea^eraturs 
Ambient pressure 
Steam to carbon ratio 
«1 Input fuel is methane 
■2 I^put fuel is methanol 
■3 loput fuel is naphtha 
ai Adiabatic operation 
m2 Isothermal operation 
Humber of coo^onents in vhole system 
Ex tra percental of needed air in burner 
Criterion of convergence in system trial and 
error procedure 

Overall beat transfer coefficient of 
reci^erator 

Transfer area of rec u perator - 

Divider fraction of input fuel 
Dsitial input fuel pressure 
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Unit 


Qefialtlon 


Vaxiable Simstislon Initial 
name value 


PPO 

6.7 

atm 

TDHDM 

394.4 

Z 

PDEUM 

2.04 

atm 

OPTC 

324 

Z 

OPPC 

4.082 

atm 

DDTYM 

0.062 


PDESIG 

1.36 

a-bs 

TCDPES 

464 

Z 

RATO 

1.8 



PDST 

0.02 

PDR 

0.1 

PDS 

0.03 

PDE 

0,005 

PDB 

0.1 

PDPO 

0.05 

PDZN 

0.03 

EFBLO 

0.75 

EFCQM 

0.6 

EEPDMP 

0.7 

X 

0.92417 

DTM2 

86.11 


Designed input fuel pressure 
Teog>eratare of steam drum 
Pressure of steam drum 
Operating ten^erature in the direct 
contact condenser 

Operating pressure in the direct contact 
condenser 

Ertra heat du-^ in steam generating 
super heater 

Designed pressure of air used in burner 
Tenperature of anode side gas reservoir 
Multiple factor of Seek corelation in 
calculating wall heat ixansfer 
coefficient in reformer 
Pressure drop ratio of steam throu^ 
steam generating super heater 
Pressure drop ratio of combustion gas 
throu^ reformer 

Presstxre drop ratio of combustion gas 
through steam generating super heater 
Pressure drop ratio of water throu^ 
economizer 

Pressure drop ratio of combustion gas 
throu^ burner 

Pressure drop ratio of anode side fuel 
in fuel cell stack 

Pressure drop ratio of input fuel throu^ 
ZnO bed 

Efficiency of air blower 

Efficiency of fuel eoopressor 

Efficiency of water pump 

Dxitial guess of conversion in reformer 

Initial guess of C02 to CO ratio in the 

anode side inlet 
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Oait Definition 


Variable Dimension Tni - h-iai 
name value 


DTM1 

0.68724 


TAI3 

806.16 

E 

PAI3 

5.159 

atm 

IA3 

374.5a 

E 

PA3 

1.3587 

atm 


Diitial guess of C02 to CO ratio after 
reformer 

Daitial guess of 13th stream teii5>erature 
Daitial guess of 1 3th stream pressure 
Initial guess of 5cd s'teeaun teo^erature 
Initial guess of 3rd s'teeam pressure 
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A6.2 Different Data Used in Chapter 6 


Variable 

name 


Dimension 


MS 

ZH 

D1 

D2 

D3 

MH 


X 

DTM2 

DTK1 

TAI3 

PAI3 

PA3 


1 

2 

1 

1 

1 

1 

1 

2 

3 

4 

5 


Data used in 
Chapter 6 


Unit 


4 

4 

2.286 

0.0762 

0.08535 

0.10973 

4 

4 

4 

4 

4 

0.923 

270. 

0.6886 

8O5.8 

5.17 

1.3599 


m 

m 

m 

m 


K 

atm 


Appendix 7 

Input Cata for Simulating Transient State Responses 

Variable Dimension Initial Unit Definition 
name value 


CRT 

5 X 10“"^ 


SDT 

3.6 X ^cr^ 

s 

MT 

9 



IDH 

9 

4,5.4,5i2 
4 , 1 , 1,3 

IFLG 

9 

9x0 

NFLG 

9 

9x0 

MFLG 

15 

15 X 2 


Criterion for testing reach of final 
steady state 

Standard (system) time interval 
Humber of major con^onents considered 
in transient simulation 
component (i) 
i=1 : air heater 

is2 : anode exhaust heat exchanger 
i=3 : preparator 

i=4 t heat exchanger- 1 (CH4 side) 

i=5 i heat exchanger-2 (air side) 

i=6 s hl^ temperature shift converter 

i=7 : low ten^erature shift converter 

i=8 : reformer 

i*9 : fuel cell stack 

Multiple of system time interval in 

component (i) 

Flaig for register of cos^ionent (i) 

Flag for operation of component (i) 
Flag for condition of gas stream 
outlets! 

2: Initial steady state 
1t transient state 
0: final steady state 


gas 

stream: 



tubeside 

of i=1 

J-2 

shellside 

of i=1 

3-3 

tubeside 

of i*2 

3-4 

shellside 

of i»2 

3-5 

tubeside 

of i=3 

3=6 

shellside 

of i=3 
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Definition 


1 f 1 tSt 3* 5* 
2.5.3.I.8. 
1 . 1 . 3 , 1.1 
9 X 1 

9x0 

15 X 0 


3=7 tubeside of i=4 

3=8 shellside of i=4 

3*9 tubeside of i*5 

3*10 shellside of i*5 

3*11 reacting gas of i*6 

3*12 reacting gas of i*7 

3*13 reforming gas of i=8 

3*14 combustion gas of i*6 - 

3*15 anode gas of i*9 

Multiple of system time interval for 

calculating hold-up of gas stream(3) 

Counted number of calcxilations in 
component (i) 

Counted system time number of component 
(i) when previous calculation performed 
lasting time number counted for 
testing final steady state 
Criterion for testing transient 
state occiirance 

Number of system time interval for 
testing final steady state 
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SUBROUTINS 


BURN 

CAlfBO 


CDPRO 

CMASS 

CMOLE 

COMP 

COHD 

CONV 

cupao 

DATACB 

DATAIN 

DITID 

DMIX 

DMIX5 

SCXiER 

EXHEJ 

exhex 

SXrijjiti 


Appendix 8 

Definitions of Subroutines and Factions 
DESCRIPTION 

calculation of mass balance In burner 

calculations of the new tesqserature of process air and cooling 
air after one time interval at transient state in fuel cell 
stack 

calculation of current density distribution at transient state 
calculation of mass fractions of gas stream 
calculation of mole fractions of gas stream 

calculatlors of power needed and outlet teiQ)erature in cosipressor 

estimation of heat du-^ in condenaer 

V/egsten method used for algebraic convergence 

calculation of steady state current density distribution on 

cell plate 

calculations of fuel input rate, process and cooling air input 
rate, and coefficients at different cooling configurations 
input data reading and change of units 
calculations of material and energy balances in divider 
calculations of material and energy balances in mixer 
calculations of material and energy balances in 3 streams* mixer 
Euler method used to solve ordinary differential equations 
estimation of the outlet temperature in tribe side with input 
heat du-fy 

estimation of the outlet tesgrerature in shell side with input 
heat du-ty 

estimation of the energy balance in heat exchanger 


SUBROOTINE 


DESCRIPTION 


flame 

HINCT 

PONCTH 

PONCTR 


PONCTS 

5OTCT2 

GAUSS 

HEPD 

LQADCE 

PROPTH 

PROPTR 


calculation of the mayimnTn fl^e ten^erature in burner 
calculations of tine rate of change for plate temperature in 
fuel cell stack 

calculations of the time rate of change for temperatures of 
tube side, shell side, and tube vail in heat exchanger 
calculations of the time rate of change for molar fraction of 
CE^, teo^erature of reforming gas, density of reforming gas, 
and teo^erature of combustion gas, vail, and catalyst in 
reformer 

calculations of the time rate of change for molar fraction of 
CO and temperatures .of reactant and catalyst in shift converter 
calculations of the nev flov rates of fuel, process air, and 
cooling air after one time interval at transient state in fuel 
cell stack 

Gauss>^eidel iteration used to solve simultaneous linear 
equations 

calculation of pressure drop in heat exchanger 

calculations of operation conditions and flov rates of fuel and 

air at transient state 

calcvilatians of the properties and coefficients of stream in 
heat exchanger after nev tes^erature and flov rate obtained 
calculations of the properties and coefficients of stream in 
reformer after nev temperature and flov rate obtain^ 
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STJBaODTINE 


DESCRIPTION 


PROPTS 

PUMP 

PVI 

SNAE 

SORT 

STEADP 

STEADH 

STEADR 

STEADS 

STMAIN 

STTRAN 

TRANC 

TRANP 


calculations of the properties aT^rf coefficients of stream in 
shift converter after new temperature and. flow rate obtained 
calculations of power and outlet ten^erature In water 

pus9 

estimations of the current density and voltage from input power 
at given conditions 

Newton^Raphson iteration used to solve equilibrium compositions 
in reformer with given equilubrium constants 
sorting of the largest input data 

calculations of steady state ourrent densi-ty and tenperature 
distributions on cell plates in fuel cell stack 
calculations of steady state solutions ( include temperature ) 
in heat exchanger 

calculations of steady state solutions ( include conversion, 
ten^eratures , and pressiares ) in reformer 
calculations of steady state solutions ( Inclimie conversion, 
ten^eratures ) in shift converter 

calculations of the steady state solutions of PARC system — 
the procedinre control subroutine 

simulation of the transient responses in PAPC system with load 
change » the proced-ure control subroutine 

execution of the simulation of heat exchanger at transient state 
execution of the simulation of fuel cell stack at transient 
state 
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SDIRODTINE 


DESCRIPTION 


TRANE 

TRANR 

VI 


execution of the simulation of heat exchanger at teansient 
state 

execution of the simulation of reformer at transient state 
calculation of relationship between voltage and current density- 
in fuel cell stack 




HJNCTION 


DESCRIPTION 


calcTilation. of heat of demethanation reaction - C02 is product 

calculation of heat of water shift reaction 

calculation of total flow rate with conversion x in reformer 

calculation of heat transfer coefficient between gas stream and 

catalyst 

calculation of heat transfer coefficient in tube side 
calculation of tube side heat transfer coefficient in reformer 
calculation of heat transfer coefficient in annular shell side 
estimation of heat transfer coefficient of annular shell side 
in reformer 

calculation of heat capacity of gas mixture 

calculation of equilibrium constant of demethanation reaction 

- CO2 is product 

calculation of equilibriiUB constant of water shift reaction 
calculation of equilibrium constant of demethanation reaction 

- CO is product 

calculation of reaction rate in low tes^erature shift converter 
calculation of reaction rate in high temperature shift converter 


HDNGE 


fourth-order Bonge-Kutta method used to solve ordinary differ- 


ential equations 


BDNGM 


Runge-Kutta^Merson method used to solve ordinary differential 


equations 

calculation of thermal conductivity of gas mixture 


calculation of total heat transfer coefficient in reformer 
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I JTOCTION DESCRIPTION 

I 

I “ ■ 

calculation of viscosity of gas mixture 

estimation equilibrium conversion of water shift reaction in 
reformar 

i 

i 

i 

f 


VIS 

X2 
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Appendix 9 

Printout of PAPC System Transient Responses 

A 9.1 Initial Steady State Solutions 
A 9.2 Transient Responses at Specific Time 
(1 sec. , 4.5 min., 9.0 min., and 13.5 min.) 

A 9.3 Pinal Steady State Solutions (14.38 min.) 
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Appendix 10 

Printout of PAFC Stack Current Density and Temperature Distributions 
at Transient State 

A 10.1 Initial Steady State 
A 10.2 Transient Distribution at 96.77 sec. 

A 10.3 Final Steady State (197.64 sec.) 


j; 


HIE WOl^^GE IS O.SSft<)20 00 A 10.1 Initial steady IDtate 







in 






UJ 












< 






Mrf 






& 





MK 

mJ 





u 

mJ 





w 

UJ 





m 

u 





UJ 






pB 

pn 





« 






mJ 

0 





B. 

z 






«c 





mJ 







UbI 





U 






V 

< 


- 









pn 

&. 





pO 






0 






e 






PM 

(9 





<v 

X 





PO 

bM 





o» 






a*-« 

0 





X * 

0 





< 0 

1© in pn 9 

© ^ pM) Pb» m 





lOM^pne 

e 0 O' iO ^ 




UJ 

u. . . . . 

e 0 9 O' ^ 




►•pn 

eo'pOi^pn 

• • • » 




<c 0 

n »9 e e 

bM«-«e e e 





• MM mm tv ¥V 





&.0 

X 


*0 *C ^ ^ 90 


pn 0 if» pn 

MM 

bMrM 09 © 

inopiMNioe 

in 0 ^ '**» 0 


^iT (VOPOp«» 

(9 in 

OmmO «9 

MM^ 0 O' oe 

9‘<90^9^9Cf^ 

^ ^ ^ iQf«. 

«0 0 ^ OP 9» «C ^ 

X MM 


e 0 e O' O' 

<v rv r*‘ <n; <v 

tv rw <v IV (V 

IV IV (V (V (V IV 

PM 9 

3 e m 0* MM 



• •••«• 


mPmM 


OmmmmmmOO 

e 

e 


c • 

^ M (V 

bM 

^ 0 

o<ne«^«ei^o 

0^««0p0b4«0 

0 0 

ib 

pM 

^ ^ 0 


M*o^i«pnrg(v 

U 

«^«(V 

< MM 9 © 9 

M .M 0 « Q 

^ PM p «0 ^ 0 

M« MM ^ M« 0 9 » 


K«OCWPP* 

X IV MM e 9 te 


pnfn«m^piiV'C 

m«npn«npn(v 

lU fO 

MM » • « • 

e e e 9 9 


' 


Oe 

<mm 99 © 

UJ 


rv 




X MM MM MM e e 

^ ^ m« i>^ 0 ^ 

^ ^ ^ ^ ^ 


UJO 

0 MM M* MM MM 

X 

fS< fO — ^ rg ^ 

pb«p^o^e^ 

in p*» ^ 0^ m pn 

X rv 

X 

PM 

9^ ^ 9^ m 0^0* • 

<NJ pn ^ *n «M 0 • 

IV pn m rv M« 0 

91V 

0 (V IV 9 m 

« ^ pn e pn pb» 

»n ^ pn K» tv 0 

pn pn po in pn po 0 

pnpnmpnpnm 

►Ml© 

^pnOiD 9 

ac<v(v«Me«c 

• •>••• 



OCIO 

< . . . . 

UJ e 0 0 e 9 




DC *M 

MM^OfV 

c. 

^ ^ e e fs) n» 

0 ^ «o e ^ 

K)rnM«<r «’0 

UJ • 

•pppiOO'O' 

XmmmmmmmmO 

*v in»n *n tv ee 

« e «n (SI «o 

9^ «n .M in M« X 

0.0 


Uj 

<* m ^ ^ fc c^ 

pn m ^ ir (V 0 

^ ^ ^ ^ tv MM MM 

c 

c 

pM 

»n ***i p»^ in tv 

pn pn Ki pn pn tn in 

pn pn pn pn tn ^ 

Ui 

uj«n»©n>» 

9 ppb pn © e 


.••«»• »M 



H- to 

bM^MMlVO' 

MM cv rj — * « 9 



0 

0 


e e e 0 9 

9> rN. m ^ u; 

p«koee*Mi^u^ 

pn MM f*» MM ^ ^ 

tbi 

UJ pn 9 * >r m 


0 «n — 0 ^ pn c 

m 0 ift e (V ^ 0 

po«omopnn> 

oc 

' 0^p^«iC 

UImMpmmmmmO 

m b*; tv e < 

iT ^ wn bn pn < 

19 m iT. in pn — # X 

<in 

^ MP M» ^ * 


pn K> po po pn po p» 

pn pn pn pn pn po b- 

PM| PP1 PO PO PT PO MM 

DC ^ 

-X 

< •- 





u; O' 

UJ o- tv IV © 

w IV 9 © « mm 

e 

0 


> *0 

> 9 p^ ^ O' 

fi. pn I*; M« e 9 

r*» ^ *n f>. art > 

fN. (sj »n e (V sp 

e»imn0>rv^ 

« MM 

« • • • • 

e e 0 e 9 

^ rw e ce <r ec 

«—fs. rv« te c 

0* iniv ^ o^ f j 

• 

«MMM PN. 9 


b"’ ^ ic b"i p*^ C5 u: 

or. 10 bn pp'i M« uj 

«n ^ m pn e»* 

UJ e 

UJ in r*» r** r- 


pn pn pmn pn pn s 

K) pp> pn p^ pn pn s 

K» ppi po po pn pn 

X 

Xp-MM — MM 



• • . • . . • b» 



p— 

W ' >r. 




296 


I 

I 



o 

<1; 

CQ 


ON 


c 

o 


3 

•H 

u 

4J 

C/3 


o Aoa 

fv IT) o 

O C^ 

lA tT) e ^ 


I I I I 
^J rg rg N 


•p 

c 

o 

•H 

w 

§ 

ii 

CVJ 

o 


X Q 
X u*i 
X CO 

. ^ rs 
u o 

C . 
tn •- o 
Cl 
rg 


o w 

00 O 


•O U kU 
15 C 
• u. < 


X O 
M X > 


Ok ^ O «o ^ 

CO -r o »o **n 

m •« %o o kTt 
rw rg rg rg rg rg 


ir> hn rg CN o o o 
^ o e CO in 
rkk CO CO #0 o m Q 

rg rg rg eg rg rg <e 


K) o m ^ ^ CO 
« C« K) *0 9> ^ o 
CO Ok o* ^ «o . 

rg rg rg rg rg rg e 


^ m •>« CO rg rg 

«o o» ^ rg pk» ^ 

Ck o o e CO <o 
rg ro ^ rg rg <n 


f> O* C* ^ ^ LJ 

rg o o Ck »n ^ O 

e ^ e o rs. < 
m ^ rg rg 


m rs o kO o %p >> 
rs rg o tn CO rs 
erg^«»C'krsuj 
m K» ro m fg fj — 


^ e ^ lO m 
to Ok ^ ro CO 
to kO to to lO kT 
rg rg rg rg rg M 


rorgo«^-g 
r^cecor^^^O 
rg rg rg rg rg rg ^ 



to 

^eeioc^^cQ 
^ «9 kp rg ^ o kO 

iO Ok ^ ^ « rs • 

rg rg rg rg rg rg e 


rg e OkT rs CO 
*n o kO rg o* PS 
^oeeror^ 
rg n K) rg rg 40 


' ro kT ^ ^ M ^ u 
^r^^Ckk^rOO 
Ok f-« ^ o Ok CO < 
rg lo fo fo rg rg ^ 


o ^ lO CO rg > 
kT o* rs wo *-N rs ■ 
o ^ o CO LJ 

*0 »0 fO fO ro pj “ 


fifiOO 
ro toco 
rg o *o tO 
to tO Ps »0 
*o<-«kr kO 
•itrgrgrg 

eeoe 

I t I I 


aooa 

«-k ^ nO «0 
o «>^ CO to 
lO kT O kT* 
kOkrpksco 
p^ 


I 1 i t 

M*«rgrg 
o e e e 

pkktfCOrOFO« OOQO 
k^OkCOCO^^ «Ps«i*jO 
lOtotOtoiOkT *««nrgo 

rgrgrgrgrgrg *-«Ok«rs 

^C0«ir *« 
eps»««ooio • • . • 
rgeokcoom oeee 

rs CO rs Pk» ps, ^ till 

rg rg rg rg rg rg 

o o o o 

to ^ fks «>« fs rg ^ 
ro^rorg-s^^OOOO 
co^^Okcopsrgrgrspso 
rgrgrgrgrgrg i ^kTkOio 

-^O^kT CO 

s>»ceaoceto 
rsio^o*o>oi»corgrOkO 
^ ^ ro rg <a o • » * • 
Okooookio^oeeo 
rgr^'Oroegfgu I i i i 



^ e e e o 
cooksOkOcookv'oeeoe 

rs. M Ok ^ o 

o> — #«oc>eo^OOOO 
rg ro «o ro rg eg o o rs rg 
• • • • » kugrspnook 

eg k^ o Ok PS 

topsseto*^rg OkO^o 
rgrsrstorgo«^k0^rg^ 
O «*<■«•<« O Ok w > • • • 
pOK>rOrOK>fj 30000 


rsOkpgkrrsCO PsOkpgkrrso 
rgrgookkTkT cvrgkO^kTkr 


orsrg^^ceiuiDr^rgOkrce 

coocopsrs«OQCcocOOrsrso 


orsrsporgio<copsrs»orgto 

egcg^pwoiOficrgcgM^or) 


rg«i«ororse&.N^OPOpso 

^okcocorspssookcocorsfs 


rs CO kT CQ rg M 
rg o> ^ o m kO 


rgoortork) 

^ Ok 10 CQ P<k rs 


•0 to rg kT rg kT 

tOfg O ^ »0 rs 


O to O rs «is to 

Ok Ok Ok CO CO rs 


rgokrgiop^ce rg^rgio^co 
eokkO*ncosOOro^sor^coo 
• ■••■•*«•••■•• 
to«ooto«o**<to*o«otoco«* 

^ O' o CO rs ps o>crocorsrs 

to 

•>cps««o 9 'Obnug«k«psp 40 >nto 

eOkkrok>>s««c.go^^Ok«>«kO 



psiooioorgarstookioerg 

OOOkcQcorst^^^OkcocOrs 


eokQMneOoo^«k«'ne 
lO ^ rs rg ^ to kT rs rg ^ ^ 

o 

C0r,^c0-^>O2c3rs-we0--r*t 

^ocrcocorsMOkCkCkceeors 


fkgcOkO'nespoeOcOkrr^Okr'ie 

Okoco^okoegokocoi^oko 


porsOkkT pskO 
tO*«^torgrg 


er co^ Ok »o PS 
0p OkOk CO CO rs 


^ 00 kO e rg to 
o rg kO r^ Pk> r*) 

e lO «•« kT 09 
e e Ok Ok 09 rs 

rgCM^M-sM 

kO CO •« *o e 
rg CO Ok kO CO r.‘ 

po ^ rg lo Ck 
e e cr Ok 09 rs 
rg rg ^ ^ ^ 

) r^ kT OO Ck o 09 
to M O kO rs (0 


p“. ^ C 9 rg Ok ^ rg ^ ^ ^ 

<^csOkcco9Psx^csOke9cOrs 


^ ^ po 09 to so es 
<ooOkOke?rsx 
«i f j rg ^ X 


PREC«B!N«- PftK BLAP^S^ MOT FSiMiD 


298 


CUULIKG OR PROCESS AIR TEMPERATURE 


r 


rw p' ^ kn o 


rj o >s 03 ^ 

9k OQ ^ fSi rs 


00 r j cj 
k."» r j 00 ^ 


«0 b*i e pxk ^ 
e^ 00 03 


P«^ fx» ^ •*? 

trtt 4 f 4 


9k 03 »0 C^ PN 
9* 9“ C^ CO 03 PS 


«0 00 ^ O 
o C“J C" ^ ^ ^ 

«o m «« 00 

O O C* C* 03 PS 

rw Pi ^ 

^ oc ^ f J ® 
rj 00 9 o 00 r« 

•o — • o rj i-** ^ 

O O Cs SA C3 PS 

rjfj — — -' — 

PO sT 00 9^ O 03 
kn •-• e «9 p« 03 

^ »0 C3 PO O 9 
O O O 9 cc *S 


m 

m 

»T 

9 >CJknpn03C' 9 'CMtf 3 *P 3039 ‘ O 

%r\ C4 m o» mi C4 inrgfO9>f0esi^ 

p*irjo3-k99'»nujp3rjo0'99*inpp><n 
9>9G0C0PspsQ£e»9>00 00rNPsOM 

•— OUJ 

eoO'9«o^pM<Oi0^iO^«4«nO 
OinP3<9Psir)QCek0P«3^Ps|0«9M 

rgv3 

oosorwco^wpsc^oooiecgooc^rs^ 

9*^9oO«OPs2I^9>9>oOOOPs^UJ 

<0 

P- eo 

r^pse^^rsps rvrs9‘^rsrs z 
e m e 9« rs o oe e u*> o o» rs e < 

on 

pp^mo^^<pp9>kno<9^p«i^ 
o es OQ ps o 9> 9> o« 00 O 

PSjMPPMPPVPAnrVMMMMPPUJ 

«n oeui 

»p*)9»%9orguj9>pn^^<0MZ« 
pppsopsrgfsiOpppsopsCMfv^Z 



P^Mfk»pvooc^p3PPpsr4koeeK< 

oe9»e»oo«oc»eo9>o»o3oougz 

rgrg-*M«ppp rgr^^pp«ppp^w 

z zc. 

p^o«p«p^oop3o^pp^ekg£ 
0 OPn« 9 MP 3 pp 00 P 3 ^ pp P 3 «p ^ Li 

O 

k9P^e3<*Ps^Z*9pnoO«9fsp«Ug 

eo99o0eop"009»^eoc30^ 

o w ^ 

^fptu*tppee-40p*ur>«poo«puip» . 
ppoeQrgp3eOL>ppooorgp3c3>Z 

uj <0 

.p“O'9C^lft03pp ^*99 bn 00^ 
PCOOC'OOOOi' O O b' C> C.3 C3 LJ U 



A, 10.3 Final Steady State 


^ KJ 00 fSI 
00 ^ o ^ ro ^ 
in ^ ^ ^ m ^ 
(V eg csj eg eg 


o tr> po oo o 

in ^ O' C' »o ^ 

m m mm ^ 

eg eg eg eg eg eg 


^ m e o 
rs.oQ 00 00 
eg eg eg eg 


egrno 
O' m 
^ m o -' 
eg eg eg ' 


^ o in 
^ pn po 
OO O' 

eg eg eg eg 


eg 

^O'OO 

o*pno 


rirje 


^ ^ ^ ^ pn po o 
eg #-4 O' O' e eg 
r< 00 r» e'* O 
eg eg eg eg eg eg eg 

eg 

m 

Of'* iO 
m lo esj ^ ^ o 
00 O' O' O' 00 r** % 
eg eg eg eg eg eg e 


^ (s: eg eg 

3 0 0 0 


cooo 

o m pn 
00 o ^ o 
pn O' « 

K> *"N po pn 

00 ^ ^ ^ 


0000 
pn m eg O' 
o m 00 in 
eg o m pn 


OOOQ 
o e o fs. 
^ ^ O CO 
eg o o 
po eg m po 
»0 O P" f" 


O 00 O P"* 00 
^ O' O' 00 pn o 
m m m m m m 
eg eg ej eg eg eg 


m pn ^ %o o pn 
^ e O' CO o 'O 
ps 00 r** r*» o 
eg eg eg eg eg eg 


ooomegego^ 
eg>rpoegeg^'OfiOOQ 
ooo' 0 ‘ 0 cop^egpop<>'oeg 
* eg eg eg eg eg eg t po oo p'* ^ 


p*. 00 00 o 
m oo eg eg 
O' e o o 
eg po fO po 


00 p'* 
eg eg*n 


*>«e m ^ ^ ^ 

eg m m eg o 00 ' - 

O'OOOO'pN.*': 

eg pn PO fo eg eg m 


0 Q^m*-^acegt^m^p'*( 

o PO po ej o ^ '• • • 

-o'oeoc'ooa^oooi 
eg 10 PO PO eg eg o I 1 1 
< 


00 e O' e 


po rO pO pO 


OQ P^ UJ - 

po^O ^ 

egegp—?'-^ 


•-1P0 o eo <0 o LU 
flO eg M O' o PO o 
O' M O O' CO < 
eg fO PO PO eg eg ^ 


•o o m «io 
o m m 
o eg ^ ^ 

pO PO pO PO 


sr' O r» % 

m ^ > Vi 
coco - 
c-* r* ui 
eg eg i: 


O' 00 o ^ ^ ’ 

eg p^ 'O m ^ ec ‘ 
o M p-« ^ e cc Lii 
10 po PO PO PO eg s: 


p'* 0 ' 0 go%er'*^moeo 4 
^ o O' rv o 
C'*-«^oo'ec*gCC£< 
egpOPOPoegcj-j'T'rrgr 

uj po PO in 4 

o c- cc m ' 

pomo'mo^ oegf"*' 

^^omej^-j'^^O'* 
o « ^ o O' igj • ♦ • 
PO PO PO PO fO eg :o o^p o 4 


i 


O' m m 'T o* 
^ o ^ oO ^ 


p" r- eg O' eo 
coco 00 ^ 


^ m P''. eg o CO 

O' O' O' o o PO 


^ompopsc 

O' O' 00 00 


00 ^ P" ^ oo 

m ^ po eg P'S 


'Opo 00 m oo *-« 

O' O' 00 oo Ps Ps 


eg ^ eg eg oo 
m m ^ 00 o 


so m o ^ e eg 

^ O' O' CO 00 Ps 


ps 0 '' ^ pO oO PO 


Ps o i-N CO ^ PO 
O' O' O' CO OC Ps 


eg eg >r ^ O' cc eg 

^ ^ gr O' €C o 

UJ 

►- O' 00 eg co ^ %e 

< O' O' O' CJ 45 Ps 


HXH COnilMG OR PROCESS AIR TEMPERATURE 



I 

t 


^ in r». iT> c* 
^ o cO ^ ^ 


r*« r J c» «o 
to eO 90 ^ rs. '<0 


m fv o 00 

^ ^ C^. O O pO 


o wn »o o 

^ O' to 90 ^ 


00 Jf ^ 90 ^ 

in ^ fM 


^ po eo aA flO ^ 
^ o^ 00 00 ^ rx. 


<M fSf rj 9» CO 
in &n ^ iO e ^ 


^ in o e ra 

< 1 ^ ^ ^ CO CO n* 


p** «>^ po CO 

^ 9 * ^ o o ^ 


^ CO ^ pO 
^ O' O' CO CO 


n/ ^ c* o eo r.* 
pM ^ o cc o 


^fM^coopo ^m^coopo 
9>r«»rsknm^ c*» r>» p*« m in ^ 


^o^poco^ouipHOOpncOpn 

9 *c^cocor«>^ 0 £C^^cOQOn^r*«> 


r^ococM 9 »pN<r^ococM 9 ^r^ 

^c^inocvirsocc^cMnocgr^ 



^«^or>^in^o«^or^pNin 

9 « 9 «c^coco^^C^^ 9 •coco^<» 


> CSf PO CM 


CMO^ o 


r^CMpnracg 


.lO^PO 
‘ ^ «0 CO 


JMCMCO 


» ^ CO 


► o *n 

> m *»r rv 


• cj in 

I C' CO 


> o CO >n 
I P** ST O 


^ *n ^ 
m 

COUJ^ 
p*n <j p-» 
•O • 
coo^^ 
^ e 
*-« csi 
or 

pnop^ 
<v ^ 
• c * 
C' s <v 
pH M o 

^ ««l f^ 

c 

000 
c^ u c^ 


pH pn ^ po pH 
^ 9» CO COPH 


(Nirumcoco 

co^«^^*n 


c» tn^^ CO 
^^c^cOf^ 


p^^opnpHtn 
pnm ^ iM 


^ ^ fu in 9 ^ 

O ^ C^C 0 P*H 
^ ^ ^ 

eocop*oo 


pHinpHpomp^ pHinpHpmn^ 
rsf^p^cocopo cM^^cocopo^ 

po^co^r^mujFO-^ooo^inpn 
9»9^cOeor^^K9^^socorHpHe 
^ ^ ^ ^ ^ ^ ^ ^ ^ 

p-> O 

fHfgOlf>PO^<pHtMOinK>^(SJ 
inoj*-<K>pHinocinfSi^»npHineM 

uJ eg 

pH^r4co^gpH&.fH.*o^sleoc4p^o 

O^^^COCOPhSC^O^C^COcOPh^ 

o 

^mooo^o ^pneco^ro 

in^cOPHn^^^ocin^copHpH^ 

• • • m • » to 

09^^0^^<00*‘^0<«^C^PN 

OCP9*9»C0rH 0^0«9^COP^ 

tn Qc 

po^%o^p4^ujK)^o^fgge3 
^rgrH»incsiC4(i^^CNipHinfgr4P» 

• • • • • -O •< 

(M*i«OCMOOC£CSI^^(M%oefi^ 

oe9»c^cocs^ee^ e^ co cc uj 

oc z. 

rVOOCJC''^*^ C'JOCCMC'KI.ir- 

c: 

srrje5Kir--<Z,r-rj«c»or^«Uj 

ooc'O'Wes — ooe>^coecC 

c cs 

^eoc^^copo^oco^^eopo^uj 
c* o c- f j eo o ^ o ^ c^ rj CO > 

1 — • 

p- in CN CO ^ in ^ C' eo ^ 
<eeCPC^coo02p' ooc^c^coeou* 

-j P 4 r>* ^ ^ • :r pg rg-^ « « ^ ~ 

C. > r- 


O' CO fj ^ ^ C* 

CP c^ O' CO CO pH 


ipHpn^o 
► ^ ^ CO 


^ PO 
pH sr o 
^ y CM 


r j pH K> ^ 

O ^ € 5 ; pH 

CJ ^ ^ ^ ^ 




301 


Tin: miriri temperature or anode side is o.a 56 A 3 d 03 k 


